V Laurain

R Tóth
email: r.toth@tudelft.nl.

W-X Zheng

M Gilson

Nonparametric Identification of LPV Models Under General Noise Conditions: An LS-SVM Based Approach

Parametric identification approaches in the Linear Parameter-Varying (LPV) setting require optimal prior selection of a set of functional dependencies, used in the parametrization of the model coefficients, to provide accurate model estimates of the underlying system. Consequently, data-driven estimation of these functional dependencies has a paramount importance, especially when very limited a priori knowledge is available. Existing overparametrization and nonparametric methods dedicated to nonlinear estimation offer interesting starting points for this problem, but need reformulation to be applied in the LPV setting. Moreover, most of these approaches are developed under quite restrictive auto-regressive noise assumptions. In this paper, a nonparametric Least-Squares Support Vector Machine (LS-SVM) approach is extended for the identification of LPV polynomial models. The efficiency of the approach in the considered noise setting is shown, but the drawback of the auto-regressive noise assumption is also demonstrated by a challenging LPV identification example. To preserve the attractive properties of the approach, but to overcome the drawbacks in the estimation of polynomial LPV models in a general noise setting, a recently developed Instrumental Variable (IV)-based extension of the LS-SVM method is applied. The performance of the introduced IV and the original LS-SVM approaches are compared in an identification study of an LPV system with unknown noise dynamics.

INTRODUCTION

Estimation of linear parameter-varying (LPV) polynomial models in an input-output (IO) setting has received a significant attention recently in the identification literature (see, e.g., [START_REF] Bamieh | Identification of linear parameter varying models[END_REF], [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], [START_REF] Hsu | Nonparametric methods for the identification of linear parameter varying systems[END_REF], [START_REF] Cerone | Set-membership identification of LPV models with uncertain measurements of the time-varying parameter[END_REF], [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]). The signal relations in these LPV models are considered to be linear just as in the linear-time invariant (LTI) case, but the parameters are assumed to be functions of a measurable time-varying signal, the so-called scheduling variable p : Z → P. Here the compact set P ⊂ R nP denotes the scheduling space. This particular structure of LPV models allows an efficient extension of the LPV prediction error framework for data-driven modeling of both nonlinear and time-varying behaviors [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]).

Most available LPV parametric identification methods, even subspace or set-membership approaches, are analyzed under the assumption that the nonlinear dependencies of the model coefficients on p are appropriately parameterized [START_REF] Bamieh | Identification of linear parameter varying models[END_REF], [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]). Such an assumption is crucial not only from the view point of stochastic consistency, but also due to issues of numerical and structural bias. The adequate selection of the p-dependent functions in a LPV model parametriza-tion is a challenging problem in general, for which the most often applied solution in the parametric case leads to a (heavy) over-parametrization of the model coefficients with a rich set of dependencies [START_REF] Tóth | Order and structural dependence selection of LPV-ARX models using a nonnegative garrote approach[END_REF]). Hence the variance of the estimates can be seriously large, even if the order of the actual model is low. This problem is similar to the classical model structure or regressor selection problem encountered in statistics and LTI identification [START_REF] Ljung | System Identification: Theory for the User[END_REF]). It also plays an important role in many nonlinear system identification problems where only a limited a priori knowledge on the model structure and hence on the involved nonlinearites is available. To cope with this problem, especially when the length of the available data-record is relatively low, an elegant solution leads through the regularization of the optimization criterion. This allows to reduce the number of estimated non-zero parameters and hence the variance on the possible expanse of bias. Following this idea, many approaches have been introduced for linear regression problems from sparse estimators to shrinkage methods. Among these methods, least-squares support vector machine's (LS-SVM's), originating form statistical learning theory [START_REF] Vapnik | Statistical Learning Theory[END_REF]), have had a significant impact on nonlinear block model identification (see [START_REF] Suykens | Least Squares Support Vector Machines[END_REF][START_REF] Falck | Linear parametric noise models for least squares support vector machines[END_REF], [START_REF] Giri | Block-oriented Nonlinear System Identification[END_REF]), due to the so-called kernel-trick, which allows a significant reduction of the variance without restricting the representation capability of the parametrization. An LS-SVM idea based identification approach has been also recently introduced in the LPV case [START_REF] Tóth | A support vector machine approach for LPV linear-regression models[END_REF]). The attractive properties of this approach have been demonstrated in terms of computational load and achieved bias and variance compared to other regularized LPV methods [START_REF] Hsu | Nonparametric methods for the identification of linear parameter varying systems[END_REF], [START_REF] Tóth | Order and structural dependence selection of LPV-ARX models using a nonnegative garrote approach[END_REF]).

However, a major drawback of regularization based approaches has been the severe restrictions on the noise, most commonly in terms of an auto-regressive noise process with the same dynamics as the deterministic system. To cope with colored noise, some modified variants of the LS-SVM algorithm has been developed in the nonlinear identification setting [START_REF] Espinoza | LS-SVM regression with autocorrelated errors[END_REF], [START_REF] Falck | Linear parametric noise models for least squares support vector machines[END_REF]), but the validity of the used assumptions on the noise and the effect of violating them have never been discussed. In practice, noise often does not satisfy completely all prior assumptions, hence bias effects can be unpredictable especially in a nonparametric context. In the parametric framework, one of the most efficient methods to handle noise modeling error is the instrumental variable (IV) technique [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]). It has been applied to many types of linear models, like LPV models [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF]), and more recently has been introduced in the LS-SVM framework [START_REF] Laurain | Introducing instrumental variables in the LS-SVM based identification framework[END_REF]).

The contribution of the current paper is to extend the LPV LS-SVM identification of polynomial IO models (see [START_REF] Tóth | A support vector machine approach for LPV linear-regression models[END_REF]) to the realistic assumption of general priori unknown noise conditions. This is achieved based on the recently introduced IV modification of the LS-SVM approach. The significance of this contribution is highlighted by an analysis of the original LS-SVM scheme under realistic (general) noise conditions, demonstrating the clear need for improvement.

The paper is organized as follows. In Section 2, the identification problem as well as the assumed model structure are introduced. The LPV LS-SVM approach is briefly summarized in Section 3. This is followed in Section 4 by an analysis of the approach via a representative example where the noise conditions violate the underlying assumption. The recently proposed IV scheme for LS-SVM methods is introduced in Section 5 and its application to LPV models is proposed in Section 6. In Section 7, both methods are compared based on their statistical properties and prediction capabilities. Finally, the conclusions are presented in Section 8.

PROBLEM DESCRIPTION

To set the preliminaries for the upcoming discussion, the concept of the data-generating system and the considered model structure are briefly introduced in this section.

The considered system

The LPV data-generating system, considered in this paper, is defined as

na i=0 a o i (p(k))y(k -i) = n b j=0 b o j (p(k))u(k -j) + v o (k), (1)
where k ∈ Z is the discrete time, u : Z → R and y : Z → R denote the input and the output signals respectively, p : Z → P is the so called scheduling variable with range P ⊆ R np , v o is a bounded stochastic noise process and, without loss of generality, a o 0 = 1. Furthermore (to keep the notation simple), the nonlinear functions a o i , b o j : P → R, describing the coefficients of the dynamic relation (1), are assumed to have static dependence on p. The latter means, that these functions depend only on p(k), i.e., the instantaneous value of p. In practice however -as it is shown in Section 7 -the results of this paper can be straightforwardly extended to the case of causal dynamic dependence (dependence on both the instantaneous and past values of p).

In the sequel, the focus of the discussion is going to be on the identification of the nonlinear functions {a o i } na i=1 and {b o j } n b j=0 , i.e., the estimation/learning of these dependencies on p. Thus, the orders n a and n b of the polynomial operators in (1) as well as the scheduling variable p are assumed to be a priori known. Moreover, both u and p are considered to be deterministic (noise-free) which is the most commonly used assumption in the LPV identification setting [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]). Furthermore, we also leave the noise v o unspecified for the time being.

The considered model

The idea upon which the intended nonparametric method is based, is to assume that the nonlinear p-dependent coefficient functions f i (p k ) (corresponding to a i and b j ) lay in a high (possibly infinite) dimensional plane:

f i (p(k)) = ng i=1 ρ ⊤ i φ i (p(k)), (2)
where each φ i : R • → R nH denotes an undefined, potentially infinite (n H = ∞) dimensional function, so-called feature map, and ρ i ∈ R nH is the associated parameter vector.

In the presented LPV context, this formulation manifests in the model structure

M ρ : y(k) = ng i=1 ρ ⊤ i φ i (p k )x i (k) + e(k), (3)
where

x i (k) = y(k -i), i = 1, . . . , n a , (4a) x na+1+j (k) = u(k -j), j = 0, . . . , n b . (4b) Additionally, introduce ρ = [ ρ ⊤ 1 . . . ρ ⊤ ng ] ⊤ ∈ R ngnH and ϕ(k) = φ ⊤ 1 (p k )x 1 (k) . . . φ ⊤ ng (p k )x ng (k) ⊤ ,
(5) such that (3) can be rewritten in the regression form

y(k) = ρ ⊤ ϕ(k) + e(k). (6)
The major difficulty here is to correctly model the nonlinearities and to minimize the risk of structural bias while keeping the computational time and the variance of the estimate relatively low. However, minimizing the structural bias by increasing the dimensions of the different φ i feature maps directly increases both the variance of the estimates and the required computational time. To cope with these issues, regularization based solutions have been reported in [START_REF] Tóth | A support vector machine approach for LPV linear-regression models[END_REF] and [START_REF] Tóth | Order and structural dependence selection of LPV-ARX models using a nonnegative garrote approach[END_REF]. In the sequel, we are going to concentrate on the nonparametric LS-SVM approach presented in [START_REF] Tóth | A support vector machine approach for LPV linear-regression models[END_REF].

THE LS-SVM METHOD FOR LPV MODELS

This section briefly summarizes the LPV LS-SVM which we intended to further extended to the general noise case. The identification criterion -the cost function by which the aim of the estimation is defined -is a regularized leastsquares (prediction) error criterion and in this approach is defined in the form of

J (ρ, e) = 1 2 ρ 2 ℓ2 + γ 2 e(t k ) 2 ℓ2 , (7)
where the bias-variance trade-off is controlled by the regularization hyper-parameter γ ∈ R + 0 . Regularization in (7) makes possible to handle large-scale over-parametrizaiton problems like in sparse estimators [START_REF] Tóth | Order and structural dependence selection of LPV-ARX models using a nonnegative garrote approach[END_REF]).

Introduce the output data vector Y = [ y(1) . . . y(N ) ] ⊤ and the regressor matrix Φ = [ ϕ(1) . . . ϕ(N ) ] ⊤ . Without going into details about the assumed noise dynamics, it can be shown that estimates a i ( ) and b j ( ) of the nonlinear functions a o i ( ) and b o j ( ) can be obtained by applying the following three steps [START_REF] Tóth | A support vector machine approach for LPV linear-regression models[END_REF]):

Step 1: Construct the matrix Ω = ΦΦ ⊤ ∈ R N ×N by using the following property

[Ω] j,k = ng i=1 [Ω (i) ] j,k (8) with [Ω (i) ] j,k = x i (j)φ ⊤ i (j)φ i (k)x i (k), = x i (j)K i p(j), p(k) x i (k).
Here, {K i } ng i=1 are positive definite so-called kernel functions. These kernels define the inner products of φ ⊤ i (j)φ i (k) and hence characterize the feature maps {φ i } ng i=1 . Specification of the kernel functions instead of φ i is called the kernel trick [START_REF] Vapnik | Statistical Learning Theory[END_REF], [START_REF] Schölkopf | Learning with Kernels[END_REF]), which allows the identification of the coefficient functions a i and b j without explicitly defining the feature maps involved. In other words, this allows the implicit definition of high dimensional feature maps without actually handling the high dimensional objects, resulting in low computational time of the estimation compared to the representation capabilities. A typical type of kernel is, for example, the Radial Basis Function (RBF) kernel:

K i (p(j), p(k)) = exp - ||p(j)-p(k)|| 2 2 σ 2 i . ( 9
)
Step 2: Use the matrix Ω to compute the dual optimal solution of (7) as

α = Ω + γ -1 I N -1 Y. (10) Step 3: Compute the estimated nonlinearities as âi ( ) = ρ⊤ i φ i ( ) = N k=1 αk x i (k)K i (p(k), ), (11a) bj ( ) = ρ⊤ j φ j ( ) = N k=1 αk x j (k)K j (p(k), ). (11b)
where j = n a + 1 + j and αk is the k th component of α. Note that the parameter vector ρ is never accessible in the SVM framework as only the combined estimates âi ( ) = ρi ⊤ φ i ( ) or bj ( ) = ρ⊤ j φ j ( ) can be computed using the kernel functions defined.

THE PITFALLS RELATED TO NOISE

So far, we have neither specified the noise v o in the data-generating system (1) nor revealed its effect on the outcome of the function estimates. Nevertheless, for an experienced reader it might have been already clear that considering the model (3) and the optimization criterion ( 7), e(k) (and v o (k) respectively) must be a white stochastic noise process to achieve unbiased estimates by the presented LS-SVM scheme. Such an assumption is often found quite restrictive in practice. Thus, to extend this estimation approach to a broader class of colored noise processes -white noise filtered by an LTI filter -, some modified variants of the LS-SVM algorithm have been developed in the nonlinear identification setting [START_REF] Espinoza | LS-SVM regression with autocorrelated errors[END_REF], [START_REF] Falck | Linear parametric noise models for least squares support vector machines[END_REF]). Nevertheless, all theses techniques have been developed by implicitly assuming that the true noise process lies in the model set defined. In order to show the weaknesses of this assumption in the LPV case, we consider an example in which the noise is also correlated with p, i.e., it has an LPV structure, a case which can often be found in LPV systems with process or actuator noise [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]). This means that an LTI assumption on the noise model is invalid.

Considered the following LPV data-generating system S o : where P ẙ and P vo are the average power of the noise-free output signal ẙ and the true noise process v o respectively.

y(k) = 2 i=1 a i (p(k), p(k -1))y(k -i) + 2 j=1 b j (p(k), p(k -1))u(k -i) + v o (k), (12a) v o (k) = 0.85sign p(k -1) -0.2 v o (k -1) + 1 + (0.2 + p(k))q -1 e o (k), ( 12b
Even knowing the correct model order, this estimation problem is challenging due to the dynamic dependency on p, the different nature of the nonlinearities involved, and the relatively low number of data with respect to the considerable power of the noise present. No a priori knowledge is supposed to be available concerning the nonlinearities. In order to avoid the selection problem of the kernels and all possible optimality considerations involved, the kernels used in the following estimation are all considered to be 2-dimensional RBF kernels with the same hyper parameter σ > 0:

K i (p j , p k ) = exp - ||pj-p k || 2 2 σ 2 , i = 1, 2, 3, 4 (13)
where p i = [ p(i) p(i -1) ] ⊤ . The aim of choosing the same kernel for all nonlinearities is to limit the number of hyperparameters to two: the regularization parameter γ and the width of the RBF's σ in (13) and hence the complexity of their adequate choice (or optimization). Note that such hyper-parameters represent a structural freedom in the LS-SVM estimation scheme, which can be commonly found in all regularization based approaches. One might argue that selection of these parameters corresponds to a model structure selection problem, which is indeed true, but with a seriously reduced degree of freedom and much better understood consequences of the choice (see [START_REF] Vapnik | Statistical Learning Theory[END_REF]). For example, increasing γ reduces the bias, but increases the variance while σ is related to the inter-sample distance and therefore the search for its adequate value can be restrained to a certain range.

In this example, the hyper-parameters γ and σ are optimized by maximizing the best fit rate (BFR) of the estimated model w.r.t. the validation data set where:

BFR = 100% • max 1 - y(k) -ŷ(k) ℓ2 y(k) -ȳ ℓ2 , 0 , (14)
with (y, u, p) being the validation data, ȳ is the mean of y and ŷ is the simulated output of the model ( 3) for (u, p), estimated by the LS-SVM on the estimation data set. This optimization has resulted in σ = 0.3 and γ = 4200 for this particular example.

In order to clearly demonstrate the pitfalls linked to the noise assumption in the LS-SVM, the "nonlinearity" in a o 1 ( ) is chosen as a simple plane. In Figure 1, the mean of the function estimates â1 ( ) provided by the LS-SVM over the 50 Monte-Carlo runs using the estimation data sets is depicted together with the true a o 1 ( ). The switching dynamical nature of the underlying noise process centered at p(k -1) = 0.2 (see (12b)) results in a clear deterioration of the function estimates, which resemble to a saturation type of nonlinearity. This bias is due to the fact that v o (k) = e o (k) which is assumed in the least-squares based estimation problem (7). Furthermore, the bias of the function estimates is different for the space p(k -1) < 0.2 and p(k -1) > 0.2. This is important because a filtering operation on the data such as proposed in [START_REF] Espinoza | LS-SVM regression with autocorrelated errors[END_REF] might suppress the bias in one of the regions, but can not efficiently suppress it on the whole scheduling domain P. This shows that in the realistic case of an unknown noise structure, even if a large number of data is available, it is not possible to separate the true underlying nature of the nonlinearity from the undesirable noise effect. This highlights a clear need for the improvement of the LPV LS-SVM method to be able to apply it under general, p-correlated noise conditions such as in the given example. Otherwise, violation of the assumed noise conditions can lead to seriously different models of the system preventing the synthesis of high-performance LPV controllers.

AN IV APPROACH FOR LS-SVM'S

To cope with general unknown noise conditions, one way is to increase the complexity of the noise model. However, this solution increases the complexity of the estimation problem, e.g., resulting in a non-convex optimization, and it does not ensure that the model is rich enough to describe the dynamics of the underlying noise process. The latter unavoidably leads to an approximation trade-off in many practical scenarios.

Another possibility is to develop an estimation method robust w.r.t. modeling error on the noise. This is the core idea behind IV approaches in the classical prediction error setting. The same concept can also be harnessed in the LS-SVM case to provide, under certain conditions, unbiased function estimates independently of the true nature of the noise [START_REF] Laurain | Introducing instrumental variables in the LS-SVM based identification framework[END_REF]).

Consider a linear-regression model which can be expressed as y(k) = ρ ⊤ ϕ(k) + e(k) and a true data-generating system which can be expressed as

y(k) = ρ ⊤ o ϕ(k) + v o (k).
The idea behind the IV approach is to introduce a so-called instrument ζ(k) ∈ R nρ in the linear-regression problem which guarantees the unbiasedness (w.r.t. the noise) of the estimate by satisfying that E{ζ(k)v o (k)} = 0 for all k ∈ Z. In other words, an IV-based method produces an unbiased estimate under the following conditions: X1 The chosen instrument ζ(k) is uncorrelated to the true noise v o (k). X2 The true noise process has a zero mean and therefore

E{v o (k)} = 0.

In the LS-SVM context, the instrument can be introduced into (6), by modifying the identification criterion (7) as

ρIV = sol 1 N N k=1 ρ + γζ(k) y(k) -ϕ ⊤ ρ = 0 . ( 15
)
It can be shown (see [START_REF] Laurain | Introducing instrumental variables in the LS-SVM based identification framework[END_REF]), that this criterion leads to unbiased function estimates by following three steps similar to the LS-SVM solution:

Step 1:

Construct the matrix Ξ = ΦZ ⊤ ∈ R N ×N where Z = [ ζ(1) . . . ζ(N ) ] ⊤ (in the LS-SVM case Ω = ΦΦ ⊤ ).
Step 2: Use the matrix Ξ to compute

α = Ξ + γ -1 I N -1 Y. ( 16
)
Step 3: The nonlinear functions f i ( ) = (ρ IV i ) ⊤ ϕ i ( ) are computed using the expression:

ρIV = N k=1 αk ζ(k) (17)
A remarkable property of this modified scheme is that it requires nearly the same computational load as the original LS-SVM approach.

THE IV-LS-SVM IN THE LPV CASE

The choice of an instrument

In LPV systems, most p-dependent polynomial noise structures (Box-Jenkins, Output Error, ARMAX, ARX, etc.) fulfill condition X2. Therefore, any instrument uncorrelated to the noise will directly imply unbiased estimates with respect to the noise. The instrument proposed below is motivated by empirical experience originating from the LTI context, but due to space restriction, its statistical properties are not studied here. The instrument proposed is defined as:

ζ(k) = φ ⊤ 1 (p k )ξ 1 (k) . . . φ ⊤ ng (p k )ξ ng (k) ⊤ , (18)
with ξ i (k) = ŷ(k -i), i = 1, . . . , n a , (19a) ξ na+1+j (k) = u(k -j), j = 0, . . . , n b , (19b)
where ŷ is either the noise-free output of the system or it denotes its estimate. As long as ŷ is uncorrelated to the noise, condition X1 is fulfilled and therefore the obtained estimate will be unbiased. The computation of ŷ will be discussed later.

Considering the instrument (18), Ξ is reads as

[Ξ] j,k = ng i=1 [Ξ i ] j,k (20)
with [Ξ (i) ] j,k = x i (j)φ ⊤ i (j)φ i (k)ξ i (k). Just as in the LPV LS-SVM case, a kernel function K i is used to define the inner product φ ⊤ i (j), φ i (k). The kernel trick is applied in the same way and therefore [Ξ (i) ] j,k = x i (j)K i p(j), p(k) ξ i (k).

(21) Using expressions ( 18) and ( 17), the estimated nonlinearities can be expressed as âi

( ) = (ρ IV i ) ⊤ φ i ( ) = N k=1 αk ξ i (k)K i (p(k), ), (22a) bj
( ) = (ρ IV j ) ⊤ φ j ( ) = N k=1 αk ξ j (k)K j (p(k), ). (22b)
where j = n a + 1 + j and αk is the k th component of α.

Computation of the instrument

In practice, the estimate ŷ of the noise-free output trajectory ŷ needs to be computed such that it is really uncorrelated to the noise v o (k). In case several realizations y i of the output are available under the same excitation conditions, it is possible to use y i as the output and y j as ŷi in the instrument for example. In the next example, a refined scheme such as in [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF] is chosen where ŷ is computed using an auxiliary model which is iteratively refined. This choice is empirically justified here and can lead to an optimal parametric estimation of LPV models as shown in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF]. The algorithm used for the Refined IV-LS-SVM method can be expressed as in Algorithm 1:

Algorithm 1 Refined IV-LS-SVM 1: Initialization: Set τ = 0 and ŷ = y. 2: repeat 3:

Compute Ξ (τ ) as in ( 21) and α(τ) using ( 16).

4:

Compute the estimated nonlinear functions â(τ) i ( ) and b(τ) j ( ) using (22a-b).

5:

Simulate the output of the obtained model ŷ(τ) and use it to generate the instrument. Increment τ by 1. 6: until ||α (τ )α (τ -1) || 2 is under a certain threshold.

COMPARISON OF THE LS AND THE IV BASED APPROACHES

To show the benefit of the proposed refined IV-LS-SVM approach in the LPV setting, the example of Section 4 is reconsidered. Using the same data sets and the same value of the hyper-parameters γ = 4200 and σ = 0.3, the estimates have been recomputed by the refined IV-LS-SVM in the 50 Monte-Carlo runs. Note that using the same γ and σ gives an advantage to the LS-SVM as these parameters were optimized for the best performance of that approach.

Similarly to the previous case, the mean of the function estimates â1 ( ) provided by the refined IV-LS-SVM over the 50 Monte-Carlo runs using the estimation data sets is depicted together with the true a o 1 ( ) in Figure 2. It clearly appears that according to the theoretical results, the bias w.r.t. the noise has completely disappeared in comparison to Figure 1. This result is not limited to the nonlinearity a o 1 . To show this, the statistical properties of the estimators are analyzed through the average bias B and the standard deviation S td of the coefficient functions on P. Using the notation āi ( . ) = Mean MC âi ( . ) for the mean and ȃi ( . ) = Std MC âi ( . ) for the standard deviation of the estimated functions over the Monte-Carlo runs, B and S td are defined as:

B = Mean P abs(a o i ( . ) -āi ( . )) , (23a)
S td = Mean P ȃi ( . ) .

Table 1 shows these values for both methods. It can be seen that the LS-SVM method has got the most bias for the nonlinearities a 1 and a 2 linked to y which is corrupted Accuracy of the estimated models in terms of the average BFR (BFR) on the validation data along the Monte-Carlo simulation is presented in Table 2. The IV-LS-SVM algorithm performs equivalently well compared to the LS-SVM method. Naturally, a performance drop occurs when the data set is further decreased as the impact of the increased variance becomes more dominant in such a scenario, reaching a critical data length, under which estimation by the IV scheme provides worse results than the LS-SVM scheme.

However in case the available data-record exceeds a critical length, the LS-SVM method cannot be directly applied on the whole data set. As the dimension N increases, the inversion of Ω ∈ N × N becomes numerically intractable. Therefore, a commonly used approach is to split the whole data set into several reasonable sized data sets and to use the mean estimated model as the final result. Therefore Figure 3 displays the evolution of the fitting score of the averaged model as the number of data points increases.

It can be seen that as expected, when the data length increases, the effect of the variance becomes negligible in comparison to the benefits brought by removing the bias.

CONCLUSION

In this a paper, two recently introduced LS-SVM based methods have been compared when they are applied to the identification of polynomial LPV models in a structural learning context. It has been shown using a representative example that, for general noise structures, the bias caused by the use of an LS-based optimization scheme can lead to serious estimation error which can mislead controller synthesis. In order to cope with the noise modeling error, a recently introduced IV-based LS-SVM method has been applied to the LPV identification problem. It has been shown that the IV-based scheme guarantees unbiased estimates at the cost of an increased variance. Nevertheless, it has been demonstrated that the proposed solution results in a relatively small variance increase even for short data records. Furthermore, the IV-LS-SVM scheme preserves the computational efficiency of the original LS-SVM scheme, which makes it an attractive identification tool for the structural exploration of LPV models.

  ) with P = [-0.3, 0.7], e o being a white noise process with zero-mean Gaussian distribution and a 1 (x, y) = -0.7 + 0.5x -0.3y, a 2 (x, y) = -0.3 + 0.4x + 0.3y -1.8xy -0.8x 2 + 0.5y 2 , b 1 (x, y) = 0.5 (sin(10x) + cos(4y)) , b 2 (x, y) = sinc(5 x 2 + y 2 ) -0.7. The available measurement data consists of 50 estimation data sets and 50 validation data sets, each with data records {u(k), y(k), p(k)} N k=1 of length N = 600, where u(k), p(k) are white stochastic processes with u(k) having a uniform distribution U(-1, 1) and p(k) ∈ U(-0.3, 0.7). The robustness of the proposed algorithms is analyzed under a signal-to-noise ratio (SNR) of 8dB, computed as SNR = 20 logP ẙ Pv o

Fig. 1 .

 1 Fig. 1. Mean estimate of a o 1 by the LS-SVM approach.

Fig. 2 .

 2 Fig. 2. Mean estimate of a o1 by the refined IV-LS-SVM approach.

  Fig. 3. Average BFR vs data length.

Table 1 .

 1 Mean and standard deviation of the estimated functions for both methods.

	Coefficient		a 1	a 2	b 1	b 2
	LS-SVM	B	0.1471	0.0838	0.0300	0.0433
		S td	0.0896	0.0935	0.0941	0.0979
	IV-LS-SVM	B	0.0291	0.0231	0.0302	0.0403
		S td	0.1062	0.1026	0.1006	0.1043

Table 2 .

 2 Average BFR on 600 data points.

	Method	BFR
	LS-SVM	73.76% ± 2.78%
	RIV-LS-SVM	76.48% ± 2.25%
	by the noise. Regarding b 1 and b 2 , both methods propose
	equivalent estimates. It can be further seen that, like in
	most other frameworks, the unbiasedness of the IV-based
	method is counterbalanced by an increased variance of
	the estimates. Nevertheless, the proposed refined scheme
	achieves a standard deviation close to the one obtained by
	the LS-SVM method.


Caroline Nadeau

Hongliang Ren

Alexandre Krupa
email: alexandre.krupa@inria.fr.

Pierre Dupont
email: pierre.dupont@childrens.harvard.edu

Target Tracking in 3D Ultrasound Volumes by Direct Visual Servoing

INTRODUCTION

Three dimensional ultrasound (3DUS) guided robotic beating heart surgery is an emerging computer-assisted surgical technique. Motion tracking of both the robot and the heart tissue is one of the key open problems in these procedures and the use of 3DUS imaging offers the possibility to extract, intra-operatively and in real time, the motion of one or more targets.

In this paper, we present a visual servoing approach to track the motion of a single target that can consist of either the surgical robot or cardiac tissues. To deal with the low quality of the 3DUS volumes, the proposed intensity-based approach requires no primitive extraction or image segmentation. This approach does not involve time consuming processing steps and can be applied to a wide range of tissue types and medical instruments.

Unlike previous work [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF] where the motion compensation task was realized physically by moving the probe attached to a robotic arm to compensate for abdominal organ motion, we propose here to track the motion of the target using a 3D region of interest (ROI) which is automatically moved within the 3D US volume to compensate for both target and probe motion. This is a new method compared with conventional 3DUS-based target motion tracking or imageless motion tracking as surveyed in [START_REF] Ren | Detection of curved robots using 3d ultrasound[END_REF]. In-vivo animal experiments were conducted to validate the tracking approach.

MATERIALS AND METHODS

3DUS guided robotic beating heart surgical procedure

Image volume sequences were acquired during an invivo porcine robotic surgery using a concentric tube robot [START_REF] Gosline | Percutaneous Intracardiac Beating-heart Surgery using Metal MEMS Tissue Approximation Tools[END_REF]. The 3D image acquisition system was a Philips IE33 system manufactured by Philips Medical (www.philips.com). The robot was inserted through a small incision in the jugular vein on the neck and then steered inside the right atrium. The surgeon positioned the matrix-array 3D probe epicardially on the heart to image the robotic instrument tip and the intracardiac tissue (Fig. 1). The system was used to acquire US volumes at a frame rate of 25 Hz for durations of 10 sec. These volumes were then loaded in a rendering software we have developed and the tracking was performed off line. Standard settings of the imaging parameters were used during image generation, including 50% overall gain, 50% compression rate, high density scan line spacing, 6 cm image depth, and zero dB power level.

Intensity-based visual servoing

In this procedure, the US probe is not actuated and motion tracking is performed in the acquired US volume by moving a ROI to compensate for probe and target motion. An image-based visual servoing strategy is considered for minimizing the error e(t) = s(t)s * between a current set of visual features s and a desired one s * by applying an instantaneous velocity v c to the ROI. To observe an exponential decrease of this visual error, the classical control law [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] is given by:

v c = -λ L s + (s(t) -s * )
, where λ is the gain involved in the exponential decrease of the error and L s + is the pseudo-inverse of an estimation of the interaction matrix L s that relates the variation of the visual features to the velocity v c .

To deal with the low quality of the 3D US images, we consider as visual features s the intensity values of the voxels of the 3D ROI:

s = {I 1,1,1 , ..., I u,v,w , ..., I L,M,N } ,
where L, M and N are respectively the width, the height and the depth of the ROI and where I u,v,w represents the intensity of the voxel of 3D coordinates (u, v, w) in the US volume. The interaction matrix associated to these intensity features has been modeled in [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF] and depends only on the 3D image gradient and the coordinates of the voxels in the US volume.

RESULTS

We conducted experiments to evaluate the ability of the algorithm to track cardiac tissue and to track the robot. The results of the soft tissue tracking task are displayed in Fig. 2. In the first volume provided by the 3D probe, the central image of the volume is displayed to allow the user to delineate the desired anatomic structure by a bounding box used as a basis for the 3D ROI considered in the control law. Then, during the subsequent volumes, the velocity computed by visual servoing is applied to this ROI. Compensated tissue motion is due predominantly to the cardiac cycle. To validate visually the tracking task, we display during the cardiac motion the interpolated US image going through the center of this ROI, first without compensation (c), then with compensation (e). In both cases, the corresponding difference images (d,f) with the desired US view are also computed. The efficiency of the compensation is shown on the difference image (f) where the ROI is roughly gray with no strong gradient, compared to the difference image (d). Moreover a visual error is defined as the Euclidean norm of the visual vector and is displayed on the curve (b). Experiments were also performed to evaluate intensity-based visual servoing for tracking robot motion during surgery. In this case, we defined the desired ROI to include the tip of the concentric tube robot, as it was navigated inside the right atrium of the beating porcine heart [START_REF] Gosline | Percutaneous Intracardiac Beating-heart Surgery using Metal MEMS Tissue Approximation Tools[END_REF]. Robot tip motion is automatically tracked by the proposed intensity approach during its movement. The results are presented in Fig. 3. In each of the acquired 3DUS volumes, the ROI remains centered on the tool tip, which validates the tracking task. Fig. 3(e) shows the decrease of the visual error with the tracking and Fig. 3(f) gives an estimate of the tool motion. This motion corresponds to a forward and backward translation along the direction of the tool shaft.

DISCUSSION

We have presented a new direct visual servoing method for the estimation and compensation of rigid motions using 3DUS volume sequences and illustrated its effectiveness using data from robotic beating-heart surgery. These experiments validated the approach both for tracking the rigid motion of robotic tools and also for tracking the quasi-rigid motion of heart tissue.

As a next step, we will investigate extending the method to deal with the non-rigid motions of heart tissue using a deformable grid instead of the rigid ROI.

Fig. 1 .

 1 Fig. 1. Beating-heart intracardiac surgery with 3DUS probe manually positioned on the epicardial surface.

Fig. 2 .

 2 Fig. 2. Tracking results for cardiac tissue. The desired anatomic target is delineated on the central view of the initial US volume (a). The motion compensation is validated by the decrease in visual error (b) and by the difference images (d,f) corresponding to the central image extracted from the ROI at iteration 240, respectively without (c) and with (e) compensation.

Fig. 3 .

 3 Fig. 3. Tracking results for the robot. The robot tip indicated by the surgeon in the initial US volume (a) is tracked in the successive ones (b,c,d). The tracking task allows the reduction of the visual error (e) and the estimation of the tool pose (f) (the θ u representation is considered to describe the orientation, where u = (u x u y u z ) ⊤ is a unit vector representing the rotation axis and θ is the rotation angle.)

Acknowledgment: This work was supported by the National Institutes of Health under grants R01HL073647 and R01HL087797 and by the ANR project US-Comp of the French National Research Agency.


Caroline Nadeau
email: caroline.nadeau@inria.fr

Alexandre Krupa

Pedro Moreira

Nabil Zemiti

Philippe Poignet

Jacques Gangloff

Active stabilization of ultrasound image for robotically-assisted medical procedures

INTRODUCTION

In the context of robotic assistance to medical gestures, we propose solutions to stabilize the ultrasound (US) image by actively compensating for the physiological motions of the patient. The considered applications are for instance the assistance for diagnoses or hepatic tumor biopsies where the liver and the tumors mainly undergo the respiratory motion. Other clinic applications, such as prostate cancer brachytherapy have been identified in [1] t h a t c o u l d b e n e f i t f r o m s u c h robotic image stabilization. The per-operative image provided by the US probe and the contact force applied to the probe are used to control the six degrees of freedom (dof) of the robot. To deal with the low quality of the US images, we propose to use the intensity of the image as visual features in a visual servoing control loop [2][3]. This vision control is associated with a force control to ensure a constant force applied to the probe. Finally, a predictive controller is implemented in the control law to take advantage of the repetitiveness of the physiological motions. We present first ex-vivo robotic results on animal tissues where we compensate for the 3D motions using successively 2D and 3D probes.

MATERIALS AND METHODS

2D and 3D intensity-based visual servoing

In order to avoid time-consuming and difficult US image processing step, we implement an intensity-based visual servoing to control the motions of the probe. This approach is applied first with a classical 2-5 MHz 2D probe (Sonosite, C60) with a frame rate of 25Hz and then with a 3D motorized one (Ultrasonix, 4DC7) that provides 5vol/s. With a 2D probe, the visual features vector is built with the intensity values of the pixels of a region of interest (ROI) of the US image. The interaction matrix L s that links the variation of these image features to the motion of the probe is modeled in [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF] and used in an imagebased visual servoing control law. This interaction matrix is computed from the 3D image gradient that is estimated once at the initial pose of the probe, without being updated during the tracking task. In the 3D case [3], a 3D deformable grid, containing a set of control points, is attached to the US volume and the intensity-based approach is applied to estimate the control points displacement between the successive volume acquired by the 3D probe. From the displacement of these points, a thin-plate spline model is considered to compute the deformation of the grid that corresponds to the non-rigid motion of the target. A rigid transformation is finally extracted from this motion model and used to control the robot by positionbased visual servoing.

Vision/force control

For guaranteeing that a constant force is applied by the probe to the skin during the tracking task, we implement an hybrid vision/force control based on an external control loop approach. In the meantime, it also improves the safety of the task. The force control is used to servo the normal translation motion of the probe while the five remaining dof are controlled by the vision.

Predictive controller

In the classic visual servoing control law, the controller does not take advantage of some knowledge on the disturbance model and more particularly of its periodicity. In this work, we implement in the control loop a repetitive predictive controller (R-GPC) in order to anticipate the effect of the disturbance [START_REF] Gangloff | Model predictive control for compensation of cyclic organ motions in teleoperated laparoscopic surgery[END_REF] .

RESULTS

The first ex-vivo robotic results on real tissues with the US intensity-based visual servoing are presented here. For this validation, a chicken stuffed with pig liver and kidneys and immersed in a water tank to avoid air gaps inside its body is carried by a 6 dof robot. A 3D periodic motion composed of combined rotations and translations is applied to this phantom. The tracking task is performed successively with a 2D and a 3D US probe, mounted on a six dof anthropomorphic robotic arm, which is also equipped with a force sensor. The experimental setup is presented in Fig. 1. Two optical markers are fixed on the probe and on the phantom and provide the relative pose of both elements thanks to an EasyTrack system. This relative pose is only used as a ground truth to validate the tracking task. The robotic arm is manually positioned above the phantom and the force control is applied with a desired force of 3N to put the probe in contact with the chicken surface. A desired ROI is then delineated in the US image and a small back and forth motion is automatically realized to acquire parallel images around the desired image and compute the interaction matrix involved in the control law. A 3D periodic motion along all translations and one rotation (around the probe axis) is applied to the phantom. The disturbance has a period of 8s a n d generates amplitudes of motion of 18mm a n d 8mm along horizontal axes and 10mm along the vertical axis with a rotation of the phantom of 4deg. The results of one tracking task with the 2D probe are displayed in Fig. 2. The ROI is initialized in the first US image (a) and the visual control is launched. The disturbance motion is applied at t=20s, then at t=160s the compensation is stopped. To visually validate the task, a visual error is defined as the Euclidean norm of this visual vector and is displayed on the curve (b). After one period of disturbance (at t=28s), this error is significantly reduced thanks to the predictive effect of the R-GPC. The curve (c) shows the respect of the force reference throughout the tracking task. Finally the curve (d), obtained thanks to the EasyTrack system, validate the robotic task in terms of pose since the relative pose of the probe with respect to the phantom is maintained constant during the tracking (t<160s). Experiments were also performed with the 3D probe. In this case, no initialization step is required to compute the image gradient and the visual servoing is launched (at t=60s) while a sinusoidal disturbance, of 12.5s of period, occurs. The results of one tracking task are presented in Fig. 3

DISCUSSION

We have presented the first ex-vivo results of actively compensated periodic motions with a robot controlled along five directions with an intensity-based visual servoing and along one direction with force. These results show the asset of the approach that does not require features extraction or segmentation of the image and is therefore well adapted for dealing with any kind of anatomic images. With a 2D probe, a frame rate of 25Hz is reached but an initialization step is required before applying the tracking. The 3D probe, on the other hand, can be used to avoid this step, but the volume rate of 5Hz limits the dynamics of the tracking. In both cases, the vision control has been coupled with a force control running at 1kHz to ensure the contact of the probe.

Fig. 1

 1 Fig. 1 Robotic setup (left) with the animal phantom (right).

Fig. 2

 2 Fig. 2 Ex-vivo results : compensation with a 2D probe.

Fig. 3

 3 Fig. 3 Ex-vivo results : compensation with a 3D probe.

Acknowledgement

This work was supported by the ANR project US-Comp of the French National Research Agency.


Maintaining visibility constraints during tele-echography with ultrasound visual servoing

Tao Li, Olivier Kermorgant and Alexandre Krupa

Abstract-This paper presents a multi-task control method to maintain the visibility of an anatomic element of interest while the doctor tele-operates a 2D ultrasound probe held by a medical robot. The prior task consists in automatically maintaining several visual constraints that guarantee an intersection between the ultrasound image plane and the anatomic object of interest and the second task allows the medical expert to manually apply probe motion through tele-operation. Unlike classical visual servoing technique which continually regulate the current visual features to desired values, our control approach gradually activates the regulation of one or several ultrasound visual features that go close to fixed limits in such a way to keep them in a safe domain. The main advantage of this approach is to give to the clinician the control of all the degrees of freedom of the probe to examine the patient while automatically preserving the visibility of the element of interest if required. Both simulations and experiments performed on an abdominal phantom demonstrate the efficiency of the visibility assistance task.

Index Terms-Tele-echography, visual servoing, visibility constraints, redundancy I. INTRODUCTION Among the different medical imaging modalities, such as ultrasound (US) imaging, X-ray imaging, magnetic resonance imaging (MRI) or Computerized Tomography (CT), the US imaging has many advantages on one hand for the doctors and on the other hand for the patients. Indeed the US modality has three major assets, it is cheap, real-time and noinvasive for the human body. Moreover the US transducer is portable and US imaging can be performed even if ferromagnetic instruments are present in the observed image. Based on these reasons, tele-operable portable ultrasound robotic system has been widely used by radiologists to perform remote ultrasound examinations of patients.

In the 90's, the tele-ultrasound examinations [START_REF] Chimiak | Architecture for a high-performance tele-ultrasound system[END_REF] required an ultrasound technician to be present by the side of the patient to manually manipulate the US probe under the order of an expert thanks to the use of a video conferencing equipment system. After that, the evolution was to actuate the US probe by a robot in such a way to automatically follow a pre-specified trajectory [START_REF] Degoulange | Hippocrate: an intrinsically safe robot for medical applications[END_REF]. However, the necessity to planify a probe trajectory reduces the flexibility of this examination procedure.

In the last decade, several research teams developed robotic systems for tele-echography as for example in [START_REF] Masuda | Three dimensional motion mechanism of ultrasound probe and its application for tele-echography system[END_REF], [START_REF] Vilchis | A new robot architecture for tele-echography[END_REF] and [START_REF] Nouaille | Optimization of a 4 dof tele-echography robot[END_REF]. With such systems, a remote radiologist can teleoperate an ultrasound probe held by a robot in such a way E-mail: Firstname.Name@inria.fr it follows the movement applyed by the expert on a haptic fictive probe. The robot "ESTELE" presented in [START_REF] Nouaille | Optimization of a 4 dof tele-echography robot[END_REF] is a light robot mounted on the patient body that was validated for abdominal and fetal echography from clinical trial conducted in 4 medical centers [START_REF] Arbeille | Realtime tele-operated abdominal and fetal echography in 4 medical centres, from one expert center, using a robotic arm & isdn or satellite link[END_REF]. However, as mentioned in [START_REF] Arbeille | Realtime tele-operated abdominal and fetal echography in 4 medical centres, from one expert center, using a robotic arm & isdn or satellite link[END_REF], a sonographer paramedic is needed by the side of the patient to hold the robot and to pre-position it on the patient in order to provide to the remote expert a first view of the organ of interest. The expert can then do the examination by teleoperating the probe. Unfortunately, sometimes the expert can lose the anatomical element of interest from the image and in this case the sonographer paramedic should re-position the robot in order to help the expert to retrieve the element of interest. For that reason, a robotic control mode can be useful to automatically maintain the visibility of the section of the element of interest in the US image while the expert tele-operates the probe. A first shared control scheme has been presented in [START_REF] Abolmaesumi | Image-guided control of a robot for medical ultrasound[END_REF] where the 3 degrees-of-freedom (DOF) of the robotized probe outside the US plane (1 translation, 2 rotations) were tele-operated and the 3 DOF lying inside the US plane (2 translations, 1 rotation) were controlled by ultrasound visual servoing to center automatically the image on the carotid artery.

In this paper, we propose another approach that gives more flexibility to the expert since all the DOF of the robot holding the probe, both in the US plane and out of the US plane, are shared between the tele-operation and an ultrasound visual task. This latter automatically activates before the observed organ of interest leaves the US observation plane in order to keep it always visible in the image with a minimum section area. This automatic assistance is also useful when the doctor exams an organ that slowly moves due to the physiological motion of the patient since it automatically adapts the probe position in such a way to always keep an intersection of the organ with the US image plane. The structure of this paper is as follows. In Section 2, we briefly present the context of our work and describe the robotic task we propose to perform in order to assist the radiologist. Section 3 presents the multi-task control law we propose. Results obtained from simulation are discussed in Section 4. Finally Section 5 demonstrates the efficiency of the robotic task from experiments performed with a robot actuating a probe in contact with a realistic abdominal phantom.

II. OBJECTIVE AND DESCRIPTION OF THE ROBOTIC TASK

Our work is realized within the ANR PROSIT French national project whose goal is to develop an interactive masterslave robotic platform for medical diagnose applications with autonomous assistance modes. The "PROSIT " robot presented in [START_REF] Nouaille | Optimization of a 4 dof tele-echography robot[END_REF] which is a modified version of "ESTELE" is shown in Fig. 1a. This 4 DOF robot can move the US probe according to: 3 rotations around X, Y and Z axes of the probe frame {R p } and 1 translation along Y-axis of the probe frame. The angular velocities around X, Y and Z axes of {R p } are denoted by ω x , ω y and ω z and v y represents the translation velocity along Y-axis of {R p }. We describe herein the successive steps of the teleechography procedure aided by our visibility assistance mode:

1) During all the procedure, the robot base is kept by an assistant on the patient body and the current US image is streamed to the expert site via internet or satellite connection.

2) The expert tele-operates all the 4 DOF of the robot with a haptic fictive probe to reach a section of the anatomic object of interest. 3) Once an organ of interest is observed, the expert can decide to activate the "visibility assistance task" by clicking inside a closed organ section S (see Fig. 1b). Its 2D contour C is then automatically extracted from the US image thanks to the real-time active contour algorithm that we presented in [START_REF] Li | A robust parametric active contour based on fourier descriptors[END_REF]. A set of visual features are then computed from C in order to track in real time the evolution of the considered organ section. 4) The doctor can then freely tele-operate all the 4 DOF unless one or several visual features go close to fixed limits. If this case occurs then the system automatically prevents the loss of the organ visibility by acting like an active virtual guide and the doctor can still teleoperate the DOF in the direction that does not degrade the visibility. 5) At any time the doctor can go back to the full teleoperation control mode to examine other organs.

III. MULTI-TASK CONTROL LAW

The classical visual servoing methods [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF] use images acquired by a camera. In this case, the 2D image is a projection of the 3D world so that the coordinates of 2D visual features are usually used to control the 6 DOF of a robotic system. However, a 2D US probe only provides information in its observed plane. Therefore, suitable features need to be used to control not only the in-plane velocity (ω z ) but also the out-of-plane velocity (ω x and ω y ) of our 4 DOF robotized probe to maintain the indicated section in sight.

We denote v = (v y , ω) the probe velocity expressed in the probe control frame, where v y is the in-plane linear velocity and ω is the angular velocity. As it is usual in ultrasound visual servo [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF], v y is regulated by a force control scheme in order to maintain a constant contact force between the probe and the patient. The visual servoing and the tele-operation regulate the angular motion ω. First, we present how the visibility constraint can be ensured with ultrasound visual servo. We then combine the visibility task with the teleoperation.

A. Ultrasound visibility constraints

There are several ways an object can leave an ultrasound image. As in classical images, the section of interest S can of course leave the field of view by passing the image border. Out-of-plane motions can also make the object disappear: in this case, the surface of S progressively becomes null and the visibility is lost. The set of visual features that handle visibility is thus defined by:

s = (x g , y g , d) (1)
where (x g , y g ) are the coordinates of the center of gravity of S and d is the square root of its surface. We take the square root instead of the actual surface in order to deal with homogeneous coordinates in s. Fig. 2 sketches the constraints that we predetermine to keep the section in the field of view. The constraints correspond to the conditions:

   x g ∈ [x -, x + ] y g ∈ [y -, y + ] d ≥ d - (2)
where [x -, x + ] is the visibility interval for x g , [y -, y + ] is the visibility interval for y g and d -is the minimum value for d. As the field of view is rendered by a circular sector, the value of x -and x + vary in function of y g to adapt its form.

A classical way to deal with the exposed constraints (2) is to define for each feature a safe interval that triggers the servoing [START_REF] Kermorgant | Combining ibvs and pbvs to ensure the visibility constraint[END_REF]. Similarly to (2), the visual feature safe region is defined as:

   x g ∈ [x s-, x s+ ] y g ∈ [y s-, y s+ ] d ≥ d s- (3)
An example of safe region is represented in Fig. 2. We now expose the control law that ensures that the defined constraints are respected.

B. A control law that copes with constraints

In classical visual servoing schemes, desired values s * are defined for each visual feature. The corresponding error e = ss * is then regulated to 0 with the following control law:

ω = -λ L + ω e ( 4
)
where λ is the control gain and L ω is angular part of the estimated interaction matrix of s expressed in the control frame [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF]. L + ω denotes the Moore-Penrose pseudo-inverse of L ω . Such control laws are known to have nice properties in terms of stability and robustness with regards to calibration imprecisions.

In our case, we want to control each visual feature only when its value approaches a constraint. The desired values are chosen in the safe region for all features: (x * g , y * g ) is defined from the center of the image and d * = d s- corresponds to the surface when the visibility task begins. In [START_REF] Kermorgant | Combining ibvs and pbvs to ensure the visibility constraint[END_REF] a weighted error, denoted He, is defined so as to ensure the visibility constraints in classical visual servoing. H is a diagonal weighting matrix that allows continuous adding or removing of visual features in the control law [START_REF] Mansard | Continuity of varying-feature-set control laws[END_REF]. In our approach, H = Diag(h x , h y , h d ) where h x , h y and h d are the varying weights that correspond to the visual features from [START_REF] Chimiak | Architecture for a high-performance tele-ultrasound system[END_REF]. Adapting the weighting formulation from [START_REF] Kermorgant | Combining ibvs and pbvs to ensure the visibility constraint[END_REF] to our visual features yields:

h x (x g ) =      xg-x s+ x + -xg if x g > x s+ xg-x s- x --xg if x g < x s- 0 else (5) h y (y g ) =      yg-y s+ y + -yg if y g > y s+ yg-y s- y --yg if y g < y s- 0 else (6) h d (d) = d-d s- d --d if d < d s- 0 else (7)
The corresponding function is represented in Fig. 3. The weight is null in the safe interval, and increases to infinite as the feature value approaches its limits. A constraint is said to be activated as soon as its weight is non-null. In [START_REF] Kermorgant | Combining ibvs and pbvs to ensure the visibility constraint[END_REF] the control law that regulates He to 0 yields:

ω = -λ(H L ω ) + He (8)
A strong hypothesis to use this control law is that H L ω should be always full rank. This is not the case in our scheme, as H L ω can be full rank 3 if all features are in the unsafe area, but may also be null if all features are in the safe region as it implies that H = 0. As we will show in simulation, using [START_REF] Li | A robust parametric active contour based on fourier descriptors[END_REF] leads to a discontinuous control law. The framework of varying-feature-set [START_REF] Mansard | Continuity of varying-feature-set control laws[END_REF] proposes a new pseudo-inverse operator, that allows using a varying-rank weighting matrix H. However, the formulation is appropriate for weights between 0 and 1. As H is defined up to a scale factor, we define the normalized weighting matrix H n as:

H n = H if H ∞ ≤ 1 H/ H ∞ else (9)
In the sequel, we denote L L Hn ω the continuous inverse of L ω with the weighting matrix H n (see [START_REF] Mansard | Continuity of varying-feature-set control laws[END_REF] for the analytical computation). Finally, the control law that ensures the visibility constraints yields:

ω = -λL L Hn ω e (10)
This control law induces no motion when all features are in their safe region. When a feature approaches a limit, the induced motion tries to make it go back to the safe region. We now explicit the interaction matrix that is used in the control law.

C. Modeling of the interaction matrix

The global interaction matrix of the US features can be written L s = L v L ω ∈ R 6×6 where L v is the linear part and L ω is the angular part that is actually used for the proposed control. Denoting s T p the twist transformation matrix that transforms velocity skews from the probe frame to the image center frame, L s yields:

L s = L s s T p (11)
where L s denotes the interaction matrix of the US features expressed in the image center. The three visual features x g , y g and d of the object section can be computed from the image moments of a segmented US section as follows:

   x g = m10 m00 y g = m01 m00 d = √ m 00 ( 12
)
where m ij is the image moments of the current section S defined by: m ij =

S

x i y j dxdy. The time variation of s can thus be expressed as follows:

ṡ =      ẋg = m00 m 2 00 ṁ10 -m10 m 2 00 ṁ00 ẏg = m00 m 2 00 ṁ01 -m01 m 2 00 ṁ00 ḋ = 1 2 √ m00 ṁ00 (13)
The interaction matrix L ij related to the image moment m ij is given in [START_REF] Mebarki | 2-d ultrasound probe complete guidance by visual servoing using image moments[END_REF] for the case of a US image section:

L ij = [m υx m υy m υz m ωx m ωy m ωz ] (14)
with:

               m υx = -i m i-1,j m υy = -j m i,j-1 m υz = x m i,j -y m i,j m ωx = x m i,j+1 -y m i,j+1 m ωy = -x m i+1,j + y m i+1,j m ωz = i m i-1,j+1 -j m i+1,j-1 (15)
where the terms x m i,j and y m i,j are computed as follows:

x m i,j = S x i y j K y dx y m i,j = S x i y j K x dy (16)

The coefficients K y and K x necessitate the knowledge of the normal vector s ∇F of the object local surface at each contour point delimiting the observed section. It is used to express the moments variation due to the out-of-plane motions of the probe. We on-line estimate this vector like in [START_REF] Mebarki | 2-d ultrasound probe complete guidance by visual servoing using image moments[END_REF] by considering successive images with their relative 3D positions obtained from the robot odometry. From ( 15) and ( 13), we can express the interaction matrix L s of the considered visual features s, expressed in the image center:

L s =     m00 m 2 00 L 10 -m10 m 2 00 L 00 m00 m 2 00 L 01 -m01 m 2 00 L 00 1 2 √ m00 L 00     (17)
From ( 11) and ( 17), the interaction matrix L ω that is used in the control law is defined. We now expose how the visibility constraints are handled during the tele-operation.

D. Fusing visibility with the haptic probe

We denote ω m = (ω xm , ω ym , ω zm ) the angular velocities that are measured from the haptic probe that is manipulated by the expert. The goal of the control law is to ensure that the actual probe velocity ω is the nearest possible to the tele-operated value ω m under the constraint of visibility.

The redundancy framework allows such hierarchical approach for the classical pseudo-inverse control law [START_REF] Li | A robust parametric active contour based on fourier descriptors[END_REF]. In this case, the global law yields:

ω = -λ(H L ω ) + He + Pω m ( 18
)
where

P = I 3 -(H L ω ) + (H L ω
) is a projection operator on the null space of the visibility task. However, in this case some of the components of ω m are totally ignored as soon as one constraint is activated, even partially. Furthermore, the behavior is discontinuous when weights are activated or deactivated. A more general operator P L Hn has been proposed in [START_REF] Mansard | A unified approach to integrate unilateral constraints in the stack of tasks[END_REF] for the varying-feature-set framework [START_REF] Mansard | Continuity of varying-feature-set control laws[END_REF]:

P L Hn = I 3 -L L Hn ω L ω ( 19
)
Such an operator allows to partially project the components of ω m when the normalized weighting matrix H n has intermediary weights. The tele-operated values ω m are thus reduced, but not ignored, when the weights are non-null. The corresponding control law yields:

ω = -λL L Hn ω e + P L
Hn ω m (20) Such a control law has been used only in highly-redundant systems [START_REF] Mansard | A unified approach to integrate unilateral constraints in the stack of tasks[END_REF] that deal with a limited number of constraints.

In our approach, the formulation of the weights allows to ensure 3 constraints while performing a tele-operated task of dimension 3, on a probe that has only 3 DOF. In the next section, simulations allow comparing control laws ( 18) and (20).

IV. SIMULATION VALIDATION

The simulator we have developed displays a binary volume which contains a sphere (Fig. 4a) to simulate the organ of interest and generates the binary images observed by a controllable 2D virtual probe. An example of binary image is shown in Fig. 4b. Red and blue lines indicate the visibility and the safe intervals. The object center of gravity is extracted and displayed in red. We now compare control laws (18) and (20) with the same predefined velocities.

The desired value for the features has been set at the image center (x * g , y * g ) = 0m and with a minimum area corresponding to d * = 0.016m. Fig. 5a and5b represent the evolution of the features during the servoing. They are normalized in order to compare easily their position in the visibility region (corresponding to [-1, +1]) and the safe region (plotted in dotted lines). In both cases, the features stay in the visibility region. x g (red) quickly reaches the upper safe value x s+ , while d (blue) reaches several times its lower safe value d s-. y g (green) oscillates in the safe region. Oscillations appear in the discontinuous control law (Fig. 5a). The corresponding weights are represented in Fig. 5c and5d. Although the weighting function may increase to infinite, small values are enough to ensure the constraints. The activation clearly corresponds to the configuration where features x g and d are near to their limits. The weighting is very different from one scheme to another, revealing very small weights for the discontinuous control law (Fig. 5c). Finally, Fig. 5e control law (Fig. 5e). On the opposite, our scheme (Fig. 5f) follows at best the predefined velocity from Fig. 4c. Changes only appear when some features are near to their limits, that is when the visibility constraints are activated.

V. ROBOTIC EXPERIMENTS

Now we present results obtained from real robotic experiments. The experimental setup is shown in Fig. 6a. An anthropomorphic 6 DOF robot arm equipped with a force sensor holds a convex 2-5 MHz US transducer that interacts with a US abdominal phantom. We purposely constraint the motion of the robot in order to simulate the 4 DOF cinematic of the "PROSIT" robot. That means that only the velocity vector v = (v y , ω x , ω y , ω z ) can be applied to the probe. A realistic abdominal phantom simulates the patient body and the automatic segmentation of the organ of interest is realtime performed (25 frames/s) by our efficient active contour algorithm that we presented in [START_REF] Li | A robust parametric active contour based on fourier descriptors[END_REF]. In the following robotic experiments, v y is regulated by the force control scheme implemented in [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF] in order to maintain a 1N contact force between the probe and the phantom. Similar desired values are defined: (x * g , y * g ) = 0m and d * = 0.021m. The main task gain is λ = 0.4.

First we move the phantom by hand (see the video accompanying this paper). No velocity is tele-operated and the goal is to ensure visibility during the phantom motion. The robot moves accordingly to the phantom motion. (bottom): constraints are activated and the object stays in the safe region. We now detail the results when tele-operated velocities are sent to the probe. The operated velocities are represented in Fig. 7b. Without control for the visibility, the surface of interest leaves the image as seen in Fig. 7a. The trajectory of the center of gravity is displayed in yellow when no constraints are taken into account and in red when at least one constraint is activated. It is displayed in white when the feature leaves the visibility region. White parts do appear in Fig. 7a at both extremities, revealing that x g leaves the acceptable interval. Furthermore, the object area is particularly small in Fig. 7a. On the opposite, Fig. 7d shows that active control prevents the object from leaving the image: the trajectory does not pass the red lines, and the object is larger in the image. Three red parts are visible: on both sides of the image, which means x g was leaving its safe interval, and in the middle of the image, that corresponds to the activation of y g . This observation is confirmed by Fig. 7e. The evolution of x g (red) shows that it first approaches its upper limit (t = 30s), then approaches its lower limit (t = 45s). This corresponds to the extremities of the trajectory. The upper limit of y g (green) is approached many times, the feature value is indeed very near to the upper safe value. Finally, the area of the surface of interest (blue) is controlled several times when approaching its lower limit. The weights are represented in Fig. 7f. The activation of each weight correspond to the US visual feature approaching its limit. The object area (blue) is activated particularly often, denoting that the operated velocity ω m would make the object surface reduce too much. The abscissa x g is activated twice, that corresponds to the two extremities already mentioned. Let us note that a weight of 0.6 is enough to ensure the visibility constraint at t = 25s. Although it is always near its upper safe value, the ordinate y g has very small weights during the scheme. The probe actual velocity (Fig. 7c) are very similar to the tele-operated ones (Fig. 7b). This shows that the user can still tele-operate the probe correctly while the visibility assistance task is active.

VI. CONCLUSION

We presented an assistance control mode that automatically maintains the visibility of an organ during teleechography examination performed with a tele-operated robot holding a 2D US probe. The main advantage of the proposed multi-task control law is that the DOF of the robotized probe are not exclusively constraint by the visibility task but also available for the tele-operation. Experimental result shows that the expert can explore the organ of interest without caring if it will leave the US view. Future works will concern the implementation of the proposed control system on the "PROSIT" robot and clinical validation during teleechography examinations.

Fig. 1 .

 1 Fig. 1. (a) The "PROSIT" tele-operable portable 4 DOF robot actuating the US probe. (b) The section S of the organ of interest observed in the 2D US image.

Fig. 2 .

 2 Fig. 2. Visibility constraints. Visibility (plain red) and safe (dotted blue) intervals are defined for xg and yg. Minimum (plain black) and safe (dotted black) values are defined for d.

Fig. 3 .

 3 Fig. 3. Weighing function. h is null in the safe interval and increases to the border.

Fig. 4 .

 4 Fig. 4. Simulation environment. 3D volume (a), 2D image (b) and operated velocity (c). The generated image observed by the virtual probe where the section of interest is segmented by the use of an active contour (displayed in red).

Fig. 5 .

 5 Fig. 5. Simulation results for control laws (18) (left) and (20) (right). Normalized visual features (top), weight values (middle) and actual probe velocity (bottom).

  Fig.6b briefly represents the feature weights (top) and the US image

Fig. 6 .

 6 Fig. 6. The experimental setup is composed by a robot holding a US probe that interacts with a realistic abdominal phantom (a). During phantom motion, the constraints are activated and the visibility is ensured (b).

  (a) US image without control. xg (white) leaves the visibility area. Surface is small.

  Probe velocity ω. ω follows ωm at best under the visibility constraint. (d) US image with control. xg goes back in the safe area. Feature weights. All features are activated during this experiment.

Fig. 7 .

 7 Fig. 7. Results on the real phantom. US image without (a) and with (d) the control law. Operated velocity ωm (b) and actual probe velocity ω (c). Normalized visual features (e) and corresponding weights (f).

ACKNOWLEDGMENT

The authors acknowledge the support of the ANR project PROSIT of the French National Research Agency and the S2E2 cluster.


Radu Vranceanu
email: vranceanu@essec.fr

Damien Besancenot
email: besancenot.damien@univ-paris13.fr

The spending multiplier in a time of massive public debt : the euro area case

Keywords: Fiscal multiplier, Euro-area, Public debt, Illiquidity, the Great Recession JEL Classi…cation: E62, G01, C23

This paper argues that in Euro-area economies, where the ECB cannot bail-out …nancially distressed governments, the …scal multiplier is adversely a¤ected by the amount of public debt. A regression model on a panel of 26 EU countries over the period 1996-2011 shows that a 10 percentage point increase in the debt-to-GDP ratio is connected to a slowdown in annual growth rates of 0.28 percentage point. Furthermore, the e¤ectiveness of …scal spending is adversely a¤ected by the amount of public debt; for a debt-to-GDP ratio above 150% the impact on growth of the …scal stimulus turns negative.

Introduction

The 2007-2009 Great Recession left behind a dramatic legacy: creeping growth and high longterm unemployment rates, a devastated banking sector and huge public debts. Three years later, Euro area countries are facing the worst di¢ culties. Financial instability and poor growth are now compounding their e¤ects in the most disruptive way. With ECB target interest rate close to zero, several economists echoed by many EU political leaders argue that only a massive …scal stimulus can revigorate the weakest economies. The basic rationale behind this view builds on traditional Keynesian reasoning: low (short-term) interest rates combined with a depressed economy signal an excess of savings over investment opportunities; in this context, a higher public spending and the connected de…cit should not push up interest rates, but stimulate growth (Krugman, 2010).

The Keynesian policy framework was developed under the assumption that bonds and money are perfect substitutes; in particular, bonds were seen as the risk-free asset. The contemporary DSGE new Keynesian macroeconomics also builds on the assumption that the central bank can maintain short and long-term interest rates on a pre-committed schedule (see [START_REF] Woodford | Simple analytics of the government expenditure multiplier, Prepared for the[END_REF].

Yet these days are gone. Public debt is no longer seen as the perfect hedge against …nancial risks, in particular in Euro area countries where the central bank cannot monetize public debt.1 If the public debt is too high, small increases in this debt, triggered by …scal stimulation of the economy, might bring about large variations in risk-adjusted interest rates. The culprit is the illiquidity risk: if investors lose con…dence in a government and refuse to hold its debt, the debt service increases and pushes the debt on an unsustainable path. 2 As an example, in 2011 and 2012, Italy, one of the four largest Euro area economies, has seen yields on 10-year Treasury bonds crossing the 7% line for several times although it is running primary surpluses. Furthermore, in "high-risk" countries, risk premia on public and corporate debt are highly correlated (IMF, 2012).

In turn, higher corporate interest rates would entail the crowding-out of private investment. If the "…nancial crowding-out" e¤ect takes over the direct public spending e¤ect, the spending multiplier can be much smaller than one and even become negative.

In this paper, we aim to provide an analysis of the spending multiplier that takes into account the …nancial risk associated to large public debts. After introducing a spending multiplier that incorporates an additional "…nancial crowding-out" e¤ect, we study the relationship between public spending, growth and public debt using a panel of 26 EU countries over the period 1996-2011.

Both growth and the growth e¤ect of a …scal stimulus appear to be weaker in high-debt countries.

Furthermore, the …scal multiplier turns negative if the public debt raises above 150% of the GDP (with a lower bound of 108% at the 95% con…dence level).

Recently several empirical papers have investigated the relationship between growth and debt in the long run. For instance [START_REF] Reinhart | Growth in a time of debt[END_REF] analyze data on forty-four countries over two hundred years. The show that the median growth rate is lower by one percentage point in countries/periods where the debt-to-GDP is above 90% as compared to countries/periods where the debt-to-GDP is below 90%. [START_REF] Kumar | Public debt and growth[END_REF] analyze a panel of advanced and emerging economies over four years; they reveal a relationship between initial debt and growth: on average, a 10 percentage point increase in initial debt-to GDP is associated to a slowdown in annual real par capita GDP growth of 0.2 percentage point per year. Similar results were obtained for euro area countries; Checherita-Westphal and Rother (2012) and [START_REF] Baum | Debt and growth: New evidence for the euro area[END_REF] show that a debtto-GDP ratio above 90-100% brings about a negative impact on growth rates. There is no explicit model of for this negative relationship between growth and debt, but, as noticed by Cottarelli and Jaramillo (2012, p.9), who reviewed this empirical literature, "high debt is expected to result in lower growth because of crowding-out e¤ects on private investment".

The paper is organized as follows. We …rstly introduce a simple analysis of the spending multiplier that incorporates the …nancial crowding out e¤ect. In the second part, we use data for EU countries in the 1996-2011 period, to gauge the impact of debt on both growth rates and the e¤ectiveness of …scal policy. The last section presents the conclusion.

2 Theory: the "augmented" spending multiplier

We start with the textbook IS condition:

Y = c[Y T (Y )] + I(i) + G;
(1)

where Y stands for the national income, c is the marginal propensity to consume (constant), with c < 1, T are taxes as a function of income, with dT =dY = T Y < 1; and G denotes public spending.

The public de…cit is DEF = G T (Y ): Private investment I is is represented as a function of i, the interest rate on corporate projects, with dI=di = I i < 0. The price level is assumed to be constant and in ‡ation expectations are zero.

Let us denote by the (subjective) probability of default on the public debt and by the risk-free interest rate (constant). Without losing much explanatory power, we can assume that in the event of default the debt residual value is zero (i.e., the haircut is 100%). With risk-neutral investors, the no-trade o¤ condition implies a risk-adjusted interest rate on public debt r de…ned by:

1 + r = 1 + 1 : (2)
It turns out that the risk-adjusted interest rate on public debt is a convex, increasing function in the subjective default probability ; r = r( ): We admit that this probability is an increasing function in the debt level D and the de…cit DEF , = F (D; DEF ); with @ =@D = F 1 (; ) > 0; and @ =@DEF = F 2 (; ) > 0: We further assume that the cross derivative @ 2 =@D@DEF = F 12 (; ) > 0 : this assumption is quite plausible: if the debt is high, chances that a higher de…cit pushes the debt out of the sustainable path should increase, and vice-versa.

In Europe, many corporations depend on public orders, subsidies and state guarantees. Under imperfect information, many investors use the country risk as a proxy for corporate risk, in particular for small, non listed …rms. Furthermore, the best risk rating of companies is in general Given the ' 12 > 0 assumption, we have d =dD < 0 : the larger the debt, the bigger the impact of a …scal stimulus on interest rate, and the more powerful is the crowding-out e¤ect.

The purpose of this elementary analysis was to emphasize the …nancial crowding-out e¤ect that can prevail when the debt-to-GDP ratio gets closer to the unsustainability limit, given investors' subjective assessment of the risk of default. The unsustainability limit has not been de…ned here, but is probably much lower for countries where the central bank cannot bail out the government, as compared to countries where the central bank can monetize public debt. A more powerful analysis would include the LM relationship, and the external sector. In an open-economy context, a higher …nancial risk would contribute to depreciate the currency, and might entail a positive e¤ect on net exports.

Empirical analysis

The European Union is made up of 27 countries in 2011. The EU-27 debt-to-GDP ratio increased from 69.9% in 1996 to 82.5% in 2011, with a strong acceleration after 2007. Only six out of the 26 countries in our sample managed to reduce the debt-to-GDP ratio. The Appendix provides a summary of the data. Besides the 11 founding members who created the Euro in 1999, six other countries joined the club later on. In the light of our analysis, this is an important institutional change insofar as at that moment they become more vulnerable to the illiquidity risk. In the second semester of 2012 the debt-to-GDP ratio reached 90.0%, compared to 71.6% at the creation of the European Monetary Union in 1999.

We estimate a growth equation, using panel on 26 EU countries over the period 1996-2011; we excluded Luxemburg given its special status and small population size. 4 Data are provided by the Eurostat online database.

Denoting by GROW T H it -the annual real growth rate of country i at time t; in percentage points; DEF it -the public de…cit as a percentage of GDP of country i at time t; DEBT it -public debt as a percentage of GDP of country i at time t; and DU M EM U -a dummy variable that takes the value 1 if the country is member of the Euro area, the growth equation is:

GROW T H it = C +a 0 DEF it +a 1 DEBT it +a 2 (DEBT it DEF it )+a 3 DU M EM U it +u i + it (6)
where u i is a country speci…c residual and it is the "usual" residual. The coe¢ cient a 2 applies to an interaction term between debt and de…cit.

We estimate equation ( 6) using both random and …xed-e¤ect panel data estimators. The panel is strongly balanced; there are 405 observations (4 observations are missing for Greece (de…cit, 1996-2001), 4 observations for Malta (growth, 1996-2001) and one observation is missing for Bulgaria in 1996). The regression model provides several interesting insights, in particular if we write the …tted growth equation ( 6) as:

GROW T H it = 6:64 0:028DEBT it + 0:006 (148 DEBT it ) DEF it 2:103DU M EM U it (7)

All things equal, a 10 percentage point increase in the debt-to-GDP ratio would bring down the growth rate by 0.28 percentage point on average, a …gure in line with the …ndings by [START_REF] Kumar | Public debt and growth[END_REF].

In line with textbook Keynesian wisdom, if the debt-to-GDP is relatively low -here, below 148% (as calculated as the ratio of …tted coe¢ cient a 0 =a 2 ) -larger de…cits (and higher spending) are associated to higher growth rates in the year when they occur. [START_REF] Hall | By how much does GDP rise if the government buys more output?[END_REF] reviewed the empirical evidence for the US and concluded that, under normal circumstances, GDP rises by roughly the amount of an increase in government purchases. According our estimates, this would happen in the EU region only if the debt-to-GDP ratio would be close to zero.

Above the critical threshold of 148% debt-to-GDP, the impact of a positive de…cit on growth rates turns negative. Notice that the 95% con…dence interval for this critical debt-to-GDP is as large as [108; 201].6

Notice also that the Euro area dummy seems to have a dramatic impact on growth. In other words, EU countries outside the Euro-area tend to outperform Euro member countries.

Of course, all these estimates should be interpreted with caution since they build essentially on linear relationships between the variables. The theoretical model has emphasized that the risk premium and the interest rate might edge up very strongly if investors believe that the debt gets closer to the sustainablity threshold.

Conclusion

Whether …scal stimulus can foster growth or not, it depends to a large extent on the response of long-term interest rates and the strength of the crowding-out e¤ect on private investment.

In this paper, we argue that the interest rate response depends on investors' assessment of the government's …nancial stability. For very large public debts, a positive de…cit and additional debt might entail a substantial rise in interest rates and a very powerful crowding-out e¤ect that can o¤set the direct e¤ect of additional spending.

We back this rationale with empirical evidence on 26 European economies over the 1996-2011 period, which includes post-recession high-debt years. It turns out that a large public debt not only has an adverse e¤ect on growth, but also reduces the positive e¤ect of a given …scal stimulus.

More in detail, if the debt-to-GDP ratio exceeds a critical threshold with a central value of 148%, the spending multiplier would become negative.

Policy implications are straightforward. For low public debt EU countries, …scal policy can be an interesting tool, to be used with cautiously. However, if high-debt European economies want to stimulate growth, the most e¢ cient solutions lie on the structural reforms side.

  Table 1 presents the output of the estimation.

	random-e¤ect model. 5
		Coe¢ cients FE Coe¢ cients RE
	C	6:640	5:991
	DEF	0:959	0:684
	DEBT	0:028	0:029
	DEBT DEF	0:006	0:005
	DU M EM U	2:103	1:522
	u	1:66	0:51
		2:93	2:93
	rho	0:24	0:03
	R 2	0:22	0:23
	Table 1. Output Estimates of the Growth Equation
	(*** = signi…cant at 1%; * = signi…cant at 10%)
	A Hausman test suggests that we can rely on the …xed-e¤ect model to properly represent
	the country-level e¤ects; for the sake of comparison we also provide the output estimate of the

The no-bail out orthodoxy of the ECB was somehow relaxed in September

2012, with introduction of the Outright Monetary Transactions programme. This new initiative allows the ECB to buy in the secondary market 3-year bonds issued by of distressed governments.2 See[START_REF] Besancenot | Default on sustainable public debt: illiquidity suspect convicted[END_REF] for a dynamic model where a small and remote risk of unsustainability can trigger imediate illiquidity default on public debt.

For instance, the correlation between the daily CDS spreads on 10-year Treasury bonds and the …ve biggest industrial corporations in Spain (a "high-risk" country), in the 2010.07-2012.07 period is 0.80. If including the two main banks, the correlation rises to 0.96.

Including Luxemburg does not change the results.

The Breusch and Pagan LM test indicates that the RE e¤ect model is better than the pooled OLS estimator.

The con…dence interval of the ratio between the …tted coe¢ cients is calculated with the Stata routine indicated at www.stata.com/statalist/archive/2010-04/msg00429.html

capped by the country rating; when a rating agency downgrades a country, the average cost of capital for corporations edges up (IMF 2012, p 39). 3 Hence, we can assume that the interest rate on corporate debt and interest rates on public debt are strongly correlated, more precisely, we assume that risk-adjusted corporate interest rate can be represented as an increasing function in r; denoted i = i(r).

Given this chain of e¤ects, we can write the corporate interest rate as an increasing function in the public debt and the de…cit:

with ' 1 > 0; ' 2 > 0; ' 12 > 0; properties that result from properties of the F function.

A …scal stimulus dG > 0 brings about an impact on de…cit, dDEF = (1 T Y dY dG )dG > 0:

Since the variation in the debt is dD = DEF; a variation in de…cit has only a second order e¤ect (positive) on the debt that we will neglect.

Then, by di¤erentiating the equilibrium condition (1):

we get the "augmented" public spending multiplier as:

The properties of the multiplier are easy to study:

If we compare with the elementary multiplier m = [1 c (1 T Y )] 1 ; it can be easily checked that < m: This is the direct consequence of the "…nancial crowding-out e¤ect".

Furthermore, if m is allays larger than one, can be smaller than one (and even negative) if the crowding out e¤ect is strong enough, which in turn depends on the response of investment to i; and the response of corporate interest rate to public debt changes: < 1 , c < ( I i ' 2 ):