Laurent Bako
email: laurent.bako@mines-douai.fr

Identification of switched linear systems via sparse optimization ⋆

Keywords: switched linear systems, system identification, sparse optimization, hybrid systems

The work presented in this paper is concerned with the identification of switched linear systems from input-output data. The main challenge with this problem is that the data are available only as a mixture of observations generated by a finite set of different interacting linear subsystems so that one does not know a priori which subsystem has generated which data. To overcome this difficulty, we present here a sparse optimization approach inspired by very recent developments from the community of compressed sensing. We formally pose the problem of identifying each submodel as a combinatorial ℓ0 optimization problem. This is indeed an NP-hard problem which can interestingly, as shown by recent literature, be relaxed into a (convex) ℓ1-norm minimization problem. We present sufficient conditions for this relaxation to be exact. The whole identification procedure allows us to extract the parameter vectors (associated with the different subsystems) one after another without any prior clustering of the data according to their respective generating-submodels. Some simulation results are included to support the potentialities of the proposed method.

Introduction

We consider in this paper the problem of identifying a switched linear system from a collection of input-output data. A switched linear system corresponds to the behavior that results from the interaction (mainly switching) between a finite set of linear dynamical subsystems. Examples of such systems arise in many different engineering fields, typically genetic regulatory networks study, air traffic management, nonlinear systems control, manufacturing processes modeling, computer vision, etc. (see e.g., [START_REF] Liberzon | Switching in systems and control[END_REF][START_REF] Sun | Switched Linear Systems: Control and Design[END_REF] for more examples). Mathematically speaking, a switched linear system can be viewed as a relation (model) of the form

y(t) = θ ⊤ λt x(t) + e(t), (1)
relating a vector x(t) ∈ R n called the regressor vector and a signal y(t) ∈ R designated as the output of the system being modeled at time t. Here, λ t ∈ {1, . . . , s} is the discrete state or the discrete mode, i.e., the index of the active submodel at time t and θ λt ∈ R n is the associated parameter vector (PV). The sequence {e(t)} of errors refers to potential mismatch or noise; it is assumed to be bounded. In a general situation, the vector x(t) appearing in [START_REF] Bako | A recursive identification algorithm for switched linear/affine models[END_REF] need not be structured but when dealing with the input-output behavior of switched dynamical systems, it sometimes takes the form

x(t) = y(t -1) • • • y(t -n a ) u(t -1) ⊤ . . . u(t -n b ) ⊤ ⊤ , (2)
where u(t) ∈ R nu and y(t) ∈ R are respectively the input and output of the considered system, n a and n b are its orders. The dimension of x(t) is therefore n = n a + n b n u and the model ( 1) is designated as a Switched Auto-Regressive eXogenous (SARX) model.

The switched system identification problem

Given observations {x(t), y(t)} N t=1 generated by a switched linear model of the form [START_REF] Bako | A recursive identification algorithm for switched linear/affine models[END_REF], with x(t) defined as in [START_REF] Bako | Algebraic identification of switched MIMO ARX models[END_REF], we are interested here in estimating the parameter vectors {θ j } s j=1 . We start by recalling from [START_REF] Vidal | Observability and identifiability of jump linear systems[END_REF] that the problem of inferring a switched model such as [START_REF] Bako | A recursive identification algorithm for switched linear/affine models[END_REF] from a set of finite measurements, admits multiple solutions so that the identification problem is not well-posed. If the structural indexes n a and n b are not fixed, then one can find for example a trivial switched linear model consisting of one single submodel with large orders that fits all the finite dataset. Even if finite and fixed values are assigned to n a and n b , there are still infinitely many switched models that explain the data. For example, it can be simply verified that there is a switched linear model with s = N submodels that can reproduce the data. In order to remove the identifiability issue, we will assume in this paper that the orders n a and n b are finite, equal for all submodels and known a priori. With this setting for the structural indexes n a and n b , the SARX of interest here will be viewed as the one that, among all switched linear models consistent with the data, has a number of submodels that is as small as possible. The interested reader is referred to the paper [START_REF] Petreczky | Identifiability of discrete-time linear switched systems[END_REF] for a more complete treatment of the identifiability problem in the framework of switched linear state space models.

Prior work

During the last ten years a number of interesting results have been achieved in the field of hybrid system identification. Examples of such works include, in the case of switched linear models, the algebraic-geometric method [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF][START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF][START_REF] Vidal | Recursive identification of switched ARX systems[END_REF], the product-of-errors based method [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF]. Other methods such as the mixed integer programming approach [START_REF] Roll | Identification of piecewise affine systems via mixed-integer programming[END_REF], the bounded-error approach [START_REF] Bemporad | A bounded-error approach to piecewise affine system identification[END_REF], the bayesian learning based procedure [START_REF] Juloski | A bayesian approach to identification of hybrid systems[END_REF], the clusteringbased strategies [START_REF] Ferrari-Trecate | A clustering technique for the identification of piecewise affine systems[END_REF][START_REF] Ferrari-Trecate | Conditions of optimal classification for piecewise affine regression[END_REF][START_REF] Nakada | Identification of piecewise affine systems based on statistical clustering technique[END_REF][START_REF] Boukharouba | Identification of piecewise affine systems based on dempster-shafer theory[END_REF] apply to piecewise affine systems, i.e., particular switched linear/affine systems where the switching surfaces are the faces of a set of non-overlapping polyhedra. An excellent survey can be found in [START_REF] Paoletti | Identification of hybrid systems: A tutorial[END_REF] where most of the methods developed prior to 2007 have been summarized. Despite the clear merits of all these pioneering contributions, one can fairly observe that the subject of hybrid system identification is still open on many challenging issues such as computational complexity reduction, optimality and convergence analysis of the proposed methods. Recently, a promising idea has emerged as to what extent some results from sparse optimization based signal recovery can be applied to hybrid system identification. This idea may indeed be an expedient for tackling simultaneously both the clustering and estimation problems that are inherent to hybrid system computation. The work of Ozay et al. [START_REF] Ozay | A sparsification approach to set membership identification of a class of affine hybrid systems[END_REF] exploits successfully such an approach. The identification of the parameter vectors is formulated as the problem of recovering a sparse vector-valued sequence, the instances of which sequence are subsequently agglomerated to reach a minimum number of submodels. The work of Elhamifar and Vidal [START_REF] Elhamifar | Sparse subspace clustering[END_REF] also suggests sparse representation as a possible alternative for solving the problem of subspace clustering. More precisely, the authors of [START_REF] Elhamifar | Sparse subspace clustering[END_REF] consider the problem of estimating bases for a set of linear/affine subspaces from data lying in the union of these subspaces. A limitation of their work however is that the mixed subspaces need to be linearly independent, an assumption which is violated when dealing with the union of more than one hyperplanes.

Contributions of this paper

The contribution of the paper consists in the development of a new identification method for switched linear systems. Data vectors generated by such systems lie in the union of a finite set of linear hyperplanes. Therefore we pose the identification of a specific submodel as the problem of extracting the hyperplane that contains the largest number of data. The corresponding submodel is hence the one that, among all submodels, achieves, over the whole dataset, the sparsest vector of fitting errors. With this formulation, one submodel can be estimated directly without any prior clustering, by means of sparse optimization, i.e., the minimization of the number of nonzero components in an error vector. Since sparse optimization is in general non-convex, it is classical to consider instead a convex ℓ 1 relaxation of this problem. We then present sufficient conditions under which the ℓ 1 relaxation is guaranteed to recover exactly the solution of the initial sparse optimization problem. In the case when these conditions are not satisfied, we show that all the PVs can still be identified by slightly adapting an iterative reweighted ℓ 1 optimization technique proposed in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. In contrast to most of the existing methods for hybrid system identification, our method lends itself to a relatively easy analysis. For example, conditions for optimality even though somewhat conservative, can be derived. A number of results from the field of compressed sensing [START_REF] Candès | Highly robust error correction by convex programming[END_REF][START_REF] Candès | Error correction via linear programming[END_REF][START_REF] Donoho | For most large underdetermined systems of linear equations the minimal ℓ 1 -norm solution is also the sparsest solution[END_REF] can be insightful for this purpose.

Outline of this paper

The remainder of the paper is organized as follows. We start by presenting in Section 2 the main mathematical terminology used in the paper. We then describe in Section 3 the proposed algorithm for the identification of switched linear systems with arbitrary switchings. Section 4 contains some numerical results that confirm the potential of our method. Section 5 concludes the paper.

Mathematical preliminaries

In this preliminary section we introduce some mathematical concepts and notations that will be extensively used throughout the paper. We first introduce a notion of k-genericity index.

Definition 1 For a given data matrix X = x(1) . . . x(N ) ∈ R n×N with n ≤ N , and for any integer k verifying 0 < k ≤ rank(X), we define the kgenericity index ν k (X) of X to be the minimum integer m such that any n × m submatrix of X has rank k:

ν k (X) = min m : ∀ (t 1 , . . . , t m ) with t i = t j for i = j, rank x(t 1 ) • • • x(t m ) = k . (3)
If k > rank(X), we set by convention ν k (X) = +∞ and if k = 0 we set ν 0 (X) = 0 for all X.

For an overview on the function ν k (•), we quickly mention the following two obvious properties.

(1) If ν k (X) = k then ν p (X) = p for all p ≤ k. (2) For any k ≤ rank(X), it holds that k ≤ ν k (X) ≤ N .
Observe additionally that when the data {x(t)} N t=1 are in general position, i.e., when any subset {x(t 1 ), . . . , x(t n )} of n data vectors are linearly independent, we have ν n (X) = n. Hence, the number ν n (X) characterizes a property of richness (or linear independence) of the columns of X. For ν n (X) to be finite, we need to assume that rank(X) = n. It will be so in all the paper. Other possible indexes for measuring linear independence between the columns of a matrix are the so-called spark and mutual coherence whose formal definitions are recalled below from [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF] and [START_REF] Donoho | Optimal sparse representation in general (nonorthogonal) dictionaries via l1 minimization[END_REF].

Definition 2

The spark of a given matrix X denoted spark(X), is the smallest number σ such that there exists a set of σ columns of X that are linearly dependent. In fact we have

spark(X) = min z∈ker(X) z =0 z 0 (4)
where ker(X) refers to the kernel subspace of matrix X and z 0 stands for the number of nonzero entries in z.

If ker(X) = {0}, then spark(X) will be conventionally set to +∞.

Definition 3

The mutual coherence of a given matrix X = [x(1), . . . , x(N )] denoted µ(X), is the largest absolute value of the cosine between different columns of X. More precisely,

µ(X) = max 1≤t,k≤N t =k x(t) ⊤ x(k) x(t) 2 • x(k) 2 (5)
where, without loss of generality 1 , it is assumed that all the columns of X are nonzero. 1 One can always remove all zero columns from X.

x 1 x 1

x 2 x 2 ν2 = 3, spark = 2. ν2 = 2, spark = 3.

O O Figure 1. Illustration of the k-genericity index, the spark and the mutual coherence for the case n = 2. The data are viewed as vectors originating from the origin O. The mutual coherence corresponds to the cosine of the smallest angle between two pairs of data vectors. For example, the mutual coherence of the data represented on the left is smaller than the one of the data on the right.

Note that the n-genericity index, the spark and the mutual coherence (see Figure 1 for a graphical illustration) are all some measures of how rich (we will say also generic) the columns of the data matrix X = [x(1) • • • x(N )] are. For example, the smaller the index ν n (X) is, the richer the data contained in X are. Similarly, small mutual coherence (or large spark) indicates highly linearly independent columns. Given a data matrix X, we may, due to the extensive use of the genericity index and the spark throughout the paper, be interested in evaluating ν n (X) or spark(X). However, using directly Definitions 1 and 2 to this purpose seems computationally difficult. We present in Lemma 4 an alternative sufficient condition for checking whether ν n (X) = n or not. Lemma 5 gives a rough idea of how large is the spark by looking at the mutual coherence, which is much easier to compute.

Lemma 4

Assume that all the columns of X ∈ R n×N are nonzero. Then for any integer p satisfying p < 1 + 1 µ(X)

, it holds that ν p (X) = p. In particular, if n <

1 + 1 µ(X) then ν n (X) = n.
PROOF. Let X be defined from X by normalizing its columns with ℓ 2 norm, i.e., the columns of X are defined as x(t) = x(t)/ x(t) 2 . Such an operation does not affect the rank properties of X so that X and X have the same genericity indexes. Our method of proof follows similar arguments as in a lemma stated in [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. Let I be a subset of {1, . . . , N } with cardinality |I| = p. Denote by XI the matrix formed with the columns of X whose indexes are in I and define G I = X⊤ I XI . According to the Gershgorin disk theorem [START_REF] Horn | Matrix analysis[END_REF], if γ is an eigenvalue of G I , then there is an i ∈ {1, . . . , p} such that

|γ -g ii | ≤ p k=1 k =i |g ik | ,
where the g ik are the entries of G I . From the definition of G I , we know that g ii = x(i) ⊤ x(i) = 1 for all i and |g ik | = x(i) ⊤ x(k) ≤ µ(X) for all (i, k) verifying i = k. We can hence write, for any eigenvalue γ of G I , |γ -1| ≤ (p -1)µ(X), which is equivalent to

-(p -1)µ(X) + 1 ≤ γ ≤ (p -1)µ(X) + 1.
Therefore if -(p -1)µ(X) + 1 > 0 or equivalently, if p < 1+1/µ(X), all the eigenvalues of G I are strictly positive. In other words G I is positive definite and rank( XI ) = rank(X I ) = p. We have hence shown that any subset of p columns of X has rank p provided p < 1 + 1/µ(X). This is equivalent to the claim of the lemma.

The following lemma is a restatement of Lemma 4 and will be of interest in the next developments.

Lemma 5 ([11]) For any matrix

X ∈ R n×N it holds that 1 + 1 µ(X) ≤ spark(X). (6)
3 The proposed solution to the identification problem

In this section we present the main results of the paper. We begin with a short overview on the switched system identification problem.

y(t) -θ ⊤ i x(t) 2 . ( 8
)
Finding an optimal solution to this problem may involve searching exhaustively over the set {1, . . . , s} N of all possible discrete state paths, an ideal path for which there is a set of corresponding PVs that fits the data. However, such a procedure is computationally very hard. Therefore, a driving principle of many approaches is to first overcome the combinatorial nature of the problem. The algebraic-geometric method [START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF] eliminates it by transforming the identification problem into the one of fitting the data with one homogeneous polynomial. In the same vein, a regularized product-of-errors criterion is formed in [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF] and optimized through continuous nonlinear optimization techniques. The recursive approach discussed in [START_REF] Bako | A recursive identification algorithm for switched linear/affine models[END_REF] alternates between data assignment to submodels and parameter vector update. The approach of the present paper is based on a formulation of the switched system identification as a sparse optimization problem which can subsequently be relaxed into a convex program.

The method

Given the data {(x(t), y(t))} N t=1 generated by the system (1), we consider the error vector

φ(θ) = y -X ⊤ θ -e (9)
where

y = [y(1) • • • y(N )] ⊤ , e = [e(1) • • • e(N )] ⊤ and X = [x(1) • • • x(N )].
Let us denote with N i the number of data (x(t), y(t)) generated exclusively by the subsystem indexed by i. Then we can observe that if θ = θ i for some i ∈ {1, . . . , s}, then φ(θ) is a sparse vector, i.e., a vector where many entries are equal to zero. More precisely, φ(θ) contains exactly N i zero entries and N -N i nonzero entries.

In order to avoid any ambiguity in the definition of the number N i , we make the following formal assumption throughout all the paper. Assumption 6 There is no data pair x(t), y(t) that fits two different submodels of the SARX (1), i.e., y(t

) = θ ⊤ i x(t) + e(t) = θ ⊤ j x(t) + e(t) =⇒ i = j.
The method to be presented relies on exploiting the sparsity of the vector φ(θ) when θ ∈ R n is a PV associated with a submodel of (1). For the sake of clarity, assume for now that the noise sequence {e(t)} is identically null and introduce the notations

θi = 1 θ ⊤ i ⊤ and x(t) = y(t) -x(t) ⊤ ⊤ . (10)
Then for any time instant t, there is i ∈ {1, . . . , s} such that

y(t) -θ ⊤ i x(t) = x(t) ⊤ θi = 0.
Hence the data record {x(t)} N t=1 lie in the union of s linear hyperplanes whose normal directions are given by the parameter vectors θi , i = 1, . . . , s. Estimating these normal vectors may require to group data lying in each hyperplane and then proceed with standard linear identification techniques for each group. Instead of doing so, we will extract the parameter vectors θ i one after another, starting directly from the entire dataset. In a sense, our method can be thought of as a robust identification approach. In fact, the method can only identify one submodel at a time and so, when identifying one submodel, data from other submodels are roughly treated as outliers or gross errors to be corrected.

The vector φ(θ) defined in ( 9) can be interpreted as the projections of the data x(t) onto a given vector θ

= 1 θ ⊤ ⊤ ∈ R n+1 , φ(θ) =      x(1) ⊤ θ . . . x(N ) ⊤ θ     . ( 11
)
To determine the PV θ i that achieves the sparsest error φ(θ i ), we can in principle solve the sparse optimization problem

min θ φ(θ) 0 , (12)
where z 0 denotes the ℓ 0 norm2 of z, that is, the number of nonzero entries of z, z 0 = |{i : z i = 0}|. Trying to solve problem ( 12) is equivalent to attempting to find a homogeneous hyperplane (or a vector θ) that contains (that is orthogonal to) as many data x(t) as possible.

If all the submodels are sufficiently excited within the data {x(t)} N t=1 then, as suggested by the following lemma, the solution to problem ( 12) is a PV representing one of the constituent submodels of system (1). Lemma 7 Assume that N i ≥ sν n (X) for all i, where s is the number of submodels in [START_REF] Bako | A recursive identification algorithm for switched linear/affine models[END_REF]. If Assumption 6 holds, then arg min

θ φ(θ) 0 = θ io ( 13
)
where i o ∈ {i : N i ≥ N j ∀j = 1, . . . , s} is one of the indices of submodels that has generated the most number of data.

PROOF. We just need to show that the solution θ = 0 to (12) lies in {θ 1 , . . . , θ s }, the true set of PVs. Suppose by contradiction that this is not the case. Then θ achieves a sparser3 error vector φ(θ) than all the θ i 's.

In other words, if we let

I(θ) = t : x(t) ⊤ θ = 0 , then |I(θ)| ≥ N i ≥ sν n (X) for all i = 1, . . . , s, where |I(θ)|
is the cardinality of I(θ). Denote with n i the number of data generated by submodel i, and whose indices are contained in I(θ). If we had n i < ν n (X) for all i, then we would get |I(θ)| = s i=1 n i < sν n (X) which would clearly violate one of the Lemma's assumptions. Therefore there is an index j such that n j ≥ ν n (X), i.e., I(θ) contains at least ν n (X) indices of data vectors x(t) generated by the same submodel j. Such data vectors form a matrix X j = x(t j 1 ) . . . x(t j nj ) ⊤ in R nj ×(n+1) . It follows that both θ and θj lie in the null space of X j , which is one dimensional by evoking the fact that n j ≥ ν n (X).

As a consequence, θ can be written as θ = λ θj with λ a nonzero scalar. Since by definition, the first entries of θ and θj are equal to 1, we have necessarily λ = 1, θ = θj and hence θ = θ j , which contradicts the initial thesis.

While Lemma 7 says that the set of parameter vectors minimizing φ(θ) 0 is included in {θ 1 , . . . , θ s }, it does not give any further information on whether the minimizer can be unique or not. Next, we characterize the uniqueness of the minimizer of [START_REF] Elhamifar | Sparse subspace clustering[END_REF] in terms of the ngenericity index of the data matrix X.

Theorem 8 If there is a vector θ satisfying

φ(θ) 0 ≤ N -ν n (X) 2 , ( 14
)
then θ is necessarily the unique vector that achieves the sparsest possible error φ(θ). If in addition, Assumption 6 holds and N ≥ (2s -1)ν n (X), then θ ∈ {θ 1 , . . . , θ s }.

PROOF. We proceed by contradiction. Assume that there is γ ∈ R n , γ = θ, achieving an error φ(γ) that is at least as sparse as φ(θ), i.e.,

φ(γ) 0 ≤ φ(θ) 0 ≤ N -ν n (X) 2
. By introducing the notation I(θ) = t : x(t) ⊤ θ = 0 , we can write

|I(θ)| = N -φ(θ) 0 ≥ N + ν n (X) 2 |I(γ)| = N -φ(γ) 0 ≥ N + ν n (X) 2 .
Using these inequalities, we can now bound the cardinality of I(θ) ∩ I(γ) as follows

|I(θ) ∩ I(γ)| = |I(θ)| + |I(γ)| -|I(θ) ∪ I(γ)| ≥ N + ν n (X) -|I(θ) ∪ I(γ)| ≥ ν n (X) (15)
because N ≥ |I(θ) ∪ I(γ)|. Let H designate the matrix formed with the vectors x(t) ∈ R n+1 whose indices t are in I(θ) ∩ I(γ). It follows from [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF] and the definition of ν n (X) that rank(H) = n and both θ and γ lie in the nullspace of H ⊤ which is of dimension one. This, together with the fact that the first entries of θ and γ are all equal to one (see e.g., Eq. ( 10) for the definition of θ from θ), lead to θ = γ, which in turn implies that θ = γ. We thus obtain a contradiction and so, the claim of the theorem holds.

The proof of the second statement follows similar steps as the proof of Lemma 7. Using the same notations, if n i were strictly less than ν n (X) for all i = 1, . . . , s, then it would hold that

1/2(N + ν n (X)) ≤ |I(θ)| = s i=1 n i < sν n (X).
It is easy to see that this is incompatible with the assumption N ≥ (2s -1)ν n (X). There is therefore at least one j such that n j > ν s (X). We can now follow the same line of arguments as in the proof of Lemma 7 to conclude that θ is necessarily in {θ 1 , . . . , θ s }.

Relaxation of the ℓ 0 problem through the basis pursuit method

Note however that the problem ( 12) is a hard non-convex optimization problem which is NP-hard in general, see e.g., [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. As a consequence, minimizing directly the cost function in ( 12) is in general intractable. A popular alternative [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF][START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] is to consider a convex relaxation of problem [START_REF] Elhamifar | Sparse subspace clustering[END_REF] based on the ℓ 1 norm. This relaxation strategy is known as the basis pursuit (BP) method and leads to the problem min

θ φ(θ) 1 , (16)
where z 1 = N i=1 |z i | for any vector z ∈ R N . This latter problem corresponds to what is classically referred to as sparse error correction problem in [START_REF] Sharon | Computation and relaxation of conditions for equivalence between ℓ 1 and ℓ 0 minimization[END_REF] and [START_REF] Candès | Error correction via linear programming[END_REF]. Contrary to the problem [START_REF] Elhamifar | Sparse subspace clustering[END_REF], the convex problem [START_REF] Horn | Matrix analysis[END_REF] can be transformed into a classical linear program which is efficiently solvable by standard convex optimization techniques [START_REF] Boyd | Convex Optimization[END_REF].

A first natural question is, under which conditions solving the convex problem ( 16) can lead to the solution of the combinatorial problem [START_REF] Elhamifar | Sparse subspace clustering[END_REF]. Second, can uniqueness be guaranteed for that solution. In order to investigate these important questions, we will slightly reformulate the optimization problem [START_REF] Horn | Matrix analysis[END_REF]. The objective is to recast it as a more standard problem known as sparse signal representation to which some results from the literature of compressed sensing [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF][START_REF] Donoho | Optimal sparse representation in general (nonorthogonal) dictionaries via l1 minimization[END_REF] might apply. Multiply Eq. ( 9) by the orthogonal projection matrix

P X = I N -X ⊤ (XX ⊤ ) -1 X. ( 17
)
Because P X X ⊤ = 0, this yields 4 P X y = P X φ(θ).

4 Remember that the noise vector e is assumed to be zero for now.

For the purpose of the analysis to be presented, we now set z = φ(θ) and replace the problem [START_REF] Horn | Matrix analysis[END_REF] with the following constrained one,

min θ,z W X z 1 subject to P X y = P X z z = y -X ⊤ θ (18)
where we have introduced a weighting matrix W X = diag( p(1) 2 , . . . , p(N ) 2 ) with the p(i) referring to columns of P X . It is assumed here that P X has no zero column so that W X is a positive definite diagonal matrix.

Hence the problem [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF] above can be viewed as an ℓ 1 relaxation of the weighted ℓ 0 problem min θ W X φ(θ) 0 which, thanks to the nonzero scale-invariance of the ℓ 0 norm, is strictly equivalent to [START_REF] Elhamifar | Sparse subspace clustering[END_REF]. Notice also that the role of W X in ( 18) is to normalize the columns of P X .

Remark 9 Note that the (N, N )-matrix P X in ( 17) can be replaced by any full row rank matrix P ∈ R (N -n)×N spanning the orthogonal complement of the column space im(X ⊤ ) of X ⊤ . Preferably, such a matrix should be selected such that it has orthogonal rows (i.e., satisfying P P ⊤ = I). The above matrix W X needs then to be defined from the columns of P .

In the following we present sufficient conditions for the convex optimization problem [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF] to uniquely recover the minimizer of [START_REF] Elhamifar | Sparse subspace clustering[END_REF]. We finally apply these results in Theorem 14 to the special context of switched linear system identification and derive sufficient conditions for determining exactly all the PVs by means of ℓ 1 norm minimization.

Theorem 10 If there is a vector θ achieving an error φ(θ) such that

φ(θ) 0 < 1 2 1 + 1 m(X) , (19)
with

m(X) = max 1≤t,k≤N t =k |M X (t, k)| (1 -M X (t, t)) (1 -M X (k, k)) M X = X ⊤ (XX ⊤ ) -1 X, ( 20
)
then θ is the unique solution to the ℓ 1 minimization problem [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF]. In Eq. ( 20), |M X (t, k)| stands for the absolute value of the (t, k)-entry of M X .

PROOF. The theorem follows directly as a consequence of Theorem 7 in [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. To see this, note that by letting z = φ(θ), problem ( 18) is equivalent to solving

min z W X z 1 subject to P X y = P X z (21)
for z and then computing θ from the linear equation z = y-X ⊤ θ. With rank(P X ) = N -n < N , this corresponds to a problem studied in [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF]. According to Theorem 7 of that paper, the vector z can be uniquely recovered if z 0 = φ(θ) 0 < 1/2 1 + 1/µ(P X ) where µ(P X ) denotes the mutual coherence of P X , see Definition 3. We are now left with showing that m(X) = µ(P X ), which follows by straightforward calculations.

Theorem 11 If there is a vector θ obeying the condition [START_REF] Liberzon | Switching in systems and control[END_REF], then θ is the unique solution to both problems [START_REF] Elhamifar | Sparse subspace clustering[END_REF] and [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF].

For the proof of Theorem 11, we need the following lemma.

Lemma 12

The orthogonal projection matrix P X defined in [START_REF] Juloski | A bayesian approach to identification of hybrid systems[END_REF] satisfies

spark(P X ) = N -ν n (X) + 1. (22)
PROOF. From the definition (4) of the spark, we know that

spark(P X ) = min z∈ker(PX ) z =0 z 0 = min z∈im(X ⊤ ) z =0 z 0 .
The second equality follows from the fact that ker(P X ) = im(X ⊤ ). To prove [START_REF] Nakada | Identification of piecewise affine systems based on statistical clustering technique[END_REF], we will, in view of the previous equalities, just show that

N -ν n (X) + 1 = min z∈im(X ⊤ ) z =0 z 0 .
By definition of the n-genericity index ν n (X), there is a subset {t 1 , . . . , t q-1 } of {1, . . . , N }, with q = ν n (X) ≥ n, such that the rank of [x(t 1 ), . . . , x(t q-1 )] is strictly less than n. As a consequence, we can find

w * ∈ R n , w * = 0, such that x(t 1 ) ⊤ w * = • • • = x(t q-1 ) ⊤ w * = 0.
For such a w * , we must have x(t) ⊤ w * = 0 for all t / ∈ {t 1 , . . . , t q-1 }. If this were not the case, i.e., if there existed a t / ∈ {t 1 , . . . , t q-1 } such that x(t) ⊤ w * = 0, then we would have w * = 0. This follows from the definition of ν n (X) which guarantees that the matrix [x(t 1 ), . . . , x(t q-1 ) x(t)]

⊤ ∈ R q×n has full column rank (since it contains ν n (X) columns). Let us now take

z * = X ⊤ w * ∈ im(X ⊤ ). Then, z * 0 = X ⊤ w * 0 is exactly equal to N -(q -1) = N -ν n (X) + 1. It follows that N -ν n (X) + 1 ≥ min z∈im(X ⊤ ) z =0 z 0 = spark(P X ).
Now let us assume by contradiction that the last inequality is strict. Then there exists at least one z ∈ im(X ⊤ ), z = 0, such that Nν n (X) + 1 > z 0 or equivalently, Nz 0 ≥ ν n (X). In words, this means that the number of zero entries in z (which is equal to Nz 0 ) is greater than ν n (X). As z ∈ im(X ⊤ ), there is w ∈ R n , w = 0 such that z = X ⊤ w. One can therefore find some time indices k 1 , . . . , k q , with q ≥ ν n (X) verifying x(k i ) ⊤ w = 0, i = 1, . . . , q. However from Definition 1 of ν n (X), q ≥ ν n (X) implies that rank ([x(k 1 ) • • • x(k q )]) = n, which in turn leads to the conclusion w = 0, contradicting the fact that z = 0.

PROOF. [Proof of Theorem 11] From Lemma 5 we know that 1 + 1 µ(PX ) ≤ spark(P X ). We have also seen in the proof of Theorem 10 that µ(P X ) = m(X) and in Lemma 12 that spark(P X ) = Nν n (X) + 1. This leads to 1 + 1 m(X) ≤ Nν n (X) + 1. Hence if there is a vector θ obeying condition [START_REF] Liberzon | Switching in systems and control[END_REF], i.e. if there is θ verifying φ(θ) 0 < 1 2 1 + 1 m(X) , then it also holds that φ(θ) 0 < 1 2 (Nν n (X) + 1). Since φ(θ) 0 is an integer, it can then be verified that φ(θ) 0 ≤ N -νn(X) 2 whether Nν n (X) + 1 is odd or even. Now, the claim of Theorem 11 follows by direct application of Theorems 8 and 10.

Remark 13 If instead of the specific matrix P X we use an arbitrary full row rank matrix P ∈ R (N -n)×N verifying P X ⊤ = 0, then Theorems 10 and 11 still hold with the number m(X) appearing in [START_REF] Liberzon | Switching in systems and control[END_REF] replaced by µ(P ). Also the equality [START_REF] Nakada | Identification of piecewise affine systems based on statistical clustering technique[END_REF] is still true if P X is replaced by P .

Theorem 11 says that if there is a vector θ producing a sufficiently sparse error φ(θ), then ( 18) is equivalent to [START_REF] Elhamifar | Sparse subspace clustering[END_REF] in the sense that they have the same (unique) solution. Since ( 18) is convex, we can therefore recover exactly the solution to the combinatorial problem [START_REF] Elhamifar | Sparse subspace clustering[END_REF] by means of convex optimization.

Without loss of generality, we can assume in this subsection that the subsystems of system (1) are indexed in such a way that N 1 ≥ N 2 ≥ . . . ≥ N s . Define X 1 = X, and for any j = 2, . . . , s, let X j be the matrix X j-1 from which all the data vectors x(t) related to the subsystem j -1 have been deleted. This way, X 1 contains N = N 1 + . . . + N s columns, X 2 contains N -N 1 columns, X 3 contains N -N 1 -N 2 columns and so forth. With these notations, we present below an immediate corollary to Theorem 11, which is relevant to the linear switched identification problem.

Theorem 14 Consider the data matrix X ∈ R n×N generated by the SARX system [START_REF] Bako | A recursive identification algorithm for switched linear/affine models[END_REF] and assume that

N 1 > N -ϑ(X 1 ) > 0, N 2 > N -N 1 -ϑ(X 2 ) > 0, . . . N s-1 > N -N 1 -• • • -N s-2 -ϑ(X s-1 ) > 0, (23)
where ϑ(X j ) = 1/2 1 + 1/m(X j ) with m(X j ) defined as in [START_REF] Ljung | System Identification Toolbox User's Guide[END_REF] for all j = 1, . . . , s. Then all the parameter vectors {θ 1 , . . . , θ s } can be extracted one after another by solving ℓ 1 minimization problems of the form [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF].

PROOF. Consider the whole dataset

X 1 = X ∈ R n×N . Then the inequality N 1 > N -ϑ(X 1 ), is equiva- lent to φ(θ 1 ) 0 = N -N 1 < ϑ(X 1
), which by Theorem 11, implies that θ 1 solves (uniquely) both ( 12) and ( 18). The result follows by applying the same reasoning to X 2 , . . . , X s-1 . Finally, X s contains data generated only by the submodel with index s; the PV θ s of this submodel can therefore be immediately obtained by ℓ 1 or ℓ 2 minimization.

As it turns out, the main difficulty in applying the basis pursuit method to computing the submodels of a hybrid system is that it may happen that none of the submodels achieves a sufficiently sparse error vector φ(θ). In other words, the conditions of Theorem 14 may not be easy to meet on arbitrary datasets. We therefore need, in more general cases, to find a way of increasing the possibility to effectively obtain, by convex ℓ 1 optimization, the solution θ which yields the sparsest error vector φ(θ).

Improving sparsity

When the PVs θ i do not realize sufficiently sparse errors φ(θ i ), we may increase the capacity of the BP method to still recover them by instead solving a weighted variant

min θ W W X φ(θ) 1 (24)
of problem ( 16) and [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF], where W = diag (w 1 , . . . , w N ) is a weighting diagonal matrix with elements w t ≥ 0.

While the role of the weighting matrix W X > 0 defined in ( 18) is to compensate potential differences of magnitude in the columns of P X 5 , W ≥ 0 in (24) designates an adjustable weighting matrix which is intended, if it is appropriately chosen, for reinforcing the sparsitypromoting ability of the ℓ 1 norm. To see this, suppose for example that we are seeking to estimate a particular PV θ = θ i for some i = 1, . . . , s. If by some means we could know (at least approximately) the entire discrete state 5 In what follows, we will not make WX appear explicitly. sequence {λ t } N t=1 , then by setting w t = 1 when λ t = i and w t = 0 when λ t = i, θ = θ i could be recovered from [START_REF] Ozay | A sparsification approach to set membership identification of a class of affine hybrid systems[END_REF]. Hence, by appropriately adjusting the weights one can encourage the obtention of the sparsest possible error vector φ(θ). Although the discrete state sequence is not known, the weights can be iteratively approached by solving a sequence of convex ℓ norm optimization problems of the form [START_REF] Ozay | A sparsification approach to set membership identification of a class of affine hybrid systems[END_REF]. This idea is supported by the results of [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] which argue that sparsity of the solution can be enhanced through reweighting ℓ 1 optimization. For the sake of completeness we recall, with a slight adaptation to our setting, the reweighted ℓ 1 optimization technique [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] in Algorithm 1 (see below).

Algorithm 1 Reweighted ℓ 1 minimization

Inputs: Data x(t), y(t)

N t=1

Initialization: Set the initial weights as: w

(0) t = 1, t = 1, . . . , N and W (0) = diag w (0) 1 , . . . , w (0) N ; Initialize a counter, r ← 0. Repeat (1) Solve the convex problem θ (r) = arg min θ W (r) φ o (θ) 1 where φ o (θ) is an ℓ 2 -normalized version of φ(θ) de- fined as φ o (θ) = x(1) ⊤ θ x(1) 2 . . . x(N ) ⊤ θ x(N ) 2 ⊤ .
(2) Update the weights as

w (r+1) t = 1 φ o t (θ (r) ) + ε , t = 1, . . . , N (3)
r ← r + 1 Until r attains a pre-specified maximum number of iterations r max or until convergence (for example when θ (r)θ (r-1) 2 < Tol, where r > 2 and Tol is a threshold).

Return θ (r)

Dealing with noisy data

We now turn to the case when the identification data are corrupted by a moderate amount of noise {e(t)}. Then recall from Eq. ( 9) that

y = X ⊤ θ + φ(θ) + e.
If the error sequence {e(t)} were known, the previous noise-free identification method could be readily applied with ỹ = ye ∈ R N as the new output vector and X ∈ R n×N as the regressor matrix. If θ represents one submodel of the switched system (1), then

φ(θ, e) = ỹ -X ⊤ θ = y -X ⊤ θ -e (25)
should be sparse, which leads to the minimization of φ(θ, e) 1 as a possible way of recovering θ. Although {e(t)} is not measurable, one can still look for the sparsest error by minimizing φ(θ) 1 and ignoring the noise e. This is only possible up to what extent the ℓ 1 optimization can stand noise. In fact the ℓ 1 estimator can correct (be insensitive to) very gross error entries with arbitrary magnitude provided the considered error vector is sparse enough. However it is much less capable of removing the effect of non-sparse small noise e than the ℓ 2 norm for example. Therefore a quite natural approach would be to minimize the ℓ 1 norm φ(θ, e) 1 of φ(θ, e) while minimizing the ℓ 2 norm e 2 of e. There are two possible ways of implementing this idea.

A first possibility is to minimize φ(θ, e) 1 while constraining e to be bounded with respect to a certain norm. In [START_REF] Candès | Highly robust error correction by convex programming[END_REF], Candès and Randall used this approach to correct errors occurring when decoding messages transmitted over communication channels. Their idea is to estimate, under of ( 25), the error e together with the vector θ which we suppose here to represent one submodel of the switched system (1). To do so, one needs however to know a priori an upper bound η on the norm of the noise. More precisely, a somewhat tight bound η satisfying e ℓ ≤ η, with ℓ a certain norm in {2, ∞, . . .}, is required. If θ is a PV for the switched linear system, then θ may be computed from the convex program min θ,e W φ(θ, e) 1 subject to e ℓ ≤ η. [START_REF] Petreczky | Identifiability of discrete-time linear switched systems[END_REF] The solution θ to ( 26) is expected to be in {θ 1 , . . . , θ s }; the vector e ∈ R N is the variable whose estimate is expected to be equal or close to the true error e defined in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. Here, W = diag (w 1 , . . . , w N ) is a positive semidefinite and diagonal weighting matrix. For both the cases ℓ = 2 and ℓ = ∞ (see [START_REF] Candès | Highly robust error correction by convex programming[END_REF] for more comments), problem ( 26) is convex and can be efficiently solved using for example interior-point or simplex methods [START_REF] Boyd | Convex Optimization[END_REF].

A second alternative approach is to trade-off φ(θ, e) 1 and e 2 2 i.e., to consider a convex optimization problem of the form

min θ,e γ φ(θ, e) 1 + 1 2 e 2 2 . (27)
In this latter method, no a priori upper bound is required on the magnitude of the noise but a regularization term γ needs to be set. Problem ( 27) can be solved for example by Iterated Shrinkage (see e.g., [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF] for some comments on these methods).

Summary of the algorithm

We have seen in the previous subsections that by applying Algorithm 1, we can identify one of the s parameter vectors of a switched system such as (1) from the whole dataset. If the conditions of Theorem 11 are satisfied, then we know that Algorithm 1 will find (after only one iteration) a vector θ * in the set {θ 1 , . . . , θ s } such that both φ(θ * ) 0 and φ(θ * ) 1 are minimum. If these conditions are not fulfilled, Algorithm 1 may not converge towards a point in {θ 1 , . . . , θ s }. However, as argued in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] and suggested by different experiments reported in [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] and in Section 4 of this paper, the algorithm is likely to find the vector θ * that realizes the sparsest error φ(θ).

According to Lemma 7 and Theorem 8, such a point θ * is in {θ 1 , . . . , θ s } when enough rich data are available.

Without loss of generality, we will denote with θ1 , i.e. the estimate of θ 1 , the point of {θ 1 , . . . , θ s } to which the algorithm converges when it is run over all the data. Observe that θ1 can be obtained from the whole mixed data without any prior clustering. Given θ1 , we need now to estimate the rest of the PVs. However we cannot proceed this time with the whole dataset because the algorithm may still converge to the same PV θ 1 . Therefore it is preferable to remove the data generated by that submodel. The indices of such data can be determined as

I( θ1 ) = t ∈ {1, . . . , N } : x(t) ⊤ θ1 x(t) 2 • θ1 2 ≤ Thres ( 28
)
where it is assumed that Tresh ∈ [0, 1] is a tolerance threshold and θ1 = 1 θ⊤ 1 ⊤ . From the data indexed by I \ I( θ1 ), we estimate θ 2 . We can repeat this procedure until all the PVs are identified. A pseudo-code of the method is summarized in Algorithm 2.

Observe that if the number of submodels is known beforehand, it can be fed into the identification Algorithm 2. In such a case the while loop can be replaced by a for loop. Also, the stopping test need not be that strict as the condition |J| = 0. The algorithm can also be stopped once the cardinality of J falls under a given number. Before pursuing this subsection, let us emphasize a few remarks.

Remark 15 (Number of submodels) When the threshold Thres in Algorithm 2 is suitably chosen, the proposed identification scheme can automatically provide the number of submodels. However when dealing with highly noisy data, selecting a threshold may, as is well-known, be a delicate task (for more details, see the paragraph concerning the implementation issues).

Algorithm 2 Identification of all PVs

(1) Inputs: {(x(t), y(t))} N t=1

(2) Initialization: S ← ∅, J ← {1, . . . , N } (3) While |J| = 0

• Estimate a submodel by the reweighted ℓ 1 minimization method (See Algorithm 1) based on the data whose indices are contained in J • Record the identified PV: S ← S ∪ {θ} • Remove from J indices of data generated by the submodel obtained:

J ← J \ (J ∩ I(θ)),
with I(θ) defined as in Eq. ( 28). ( 4) EndWhile (5) Return S and s = |S|.

Remark 16 (Unknown orders) In case the orders n a and n b are unknown, assume that some upper bounds na ≥ n a and nb ≥ n b are available a priori. We can then parameterize the regressor x(t) with na and nb i.e., define a regressor

x na,nb (t) = y(t -1) • • • y(t -na ) u(t -1) ⊤ . . . u(t -nb ) ⊤ ⊤ .
In this case, we can, instead of [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF], solve the convex problem min

θ W φ φ(θ) 1 + W θ θ 1
where W φ and W θ are some diagonal weighting matrices with positive entries. This way, the solution θ is expected to be sparse; its sparsity can be enhanced by choosing a matrix W θ so as to penalize more the entries of θ which correspond to the components y(ti) and u(tj) of x na,nb (t) having the largest i and j.

Implementation issues. Implementation of the identification Algorithm 2 necessitates mainly three userdefined parameters to be set. The number ε in Algorithm 1 aims essentially to prevent division by zero; the tolerance Tol in the same algorithm is used to detect convergence. Finally, the number Thres appearing in Algorithm 2 is a decision parameter. As is apparent from their respective roles, the first two user-defined parameters are easy to tune and do not have much effect on the performance of the algorithm. The parameter Tol can even be replaced with a fixed maximum number of iterations (See Algorithm 1). The last parameter is probably the most delicate to set when the identification data are affected by noise. The difficulty in choosing Thres increases with the magnitude of the noise. First notice that the threshold Thres does not impact the identification of the first PV θ 1 (first is used here with respect to the numbering agreed in the beginning of this subsection) because all the dataset is processed at this first step. However, it does have an influence on the identification of the subsequent PVs. For example the identification of the second PV θ 2 requires that we manipulate the dataset in the objective of promoting the convergence of the algorithm to that second PV. An intuitive way of achieving this goal is, as proposed in Algorithm 2, to remove the data generated by the first submodel, before proceeding further. These data need not be removed entirely. For the algorithm to choose θ 2 as a convergence point, we just need a sufficient number of them to be removed. What can happen however is that if Thres is too small with respect to the level of noise, only a very small number of the data related to submodel 1 will be removed so that the algorithm can still identify θ 1 in the second step. This will likely results in an overestimation of the number of submodels with redundant PVs. Although the number of estimated submodels is not minimal, the input-output map can still be reconstructed from the identified PVs. If Thres is too large, more than the data pertaining to the first submodel may be removed, probably causing the algorithm to provide bad estimates for some of the remaining PVs.

Applications

In this section, we apply the proposed identification algorithm to a SISO SARX model composed of three linear submodels of order two. The SARX model is defined by

y(t) = θ ⊤ λt y(t -1) y(t -2) u(t -1) u(t -2) ⊤ + e(t) (29
Using this switched model, we generate the identification data under the following conditions:

• The excitation input {u(t)} is a centered signal with normal distribution and variance unity. • The noise {e(t)} is a white Gaussian noise whose magnitude is such that the Signal to Noise Ratio (SNR) is equal to 30 dB with respect to the output signal. • The switching sequence {λ t } is uniformly distributed in {1, 2, 3}.

For all the results to be presented in this section, all convex optimization instances occurring in our method have been numerically implemented with the Matlab based software developed by Grant and Boyd [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF].

Tests on the BP and the reweighted ℓ 1 methods

In a first experiment, we verify the ability of the BP method to exactly recover the solution of the sparse optimization problem, under the condition derived in Theorem 11. With regard to this goal we can set the noise sequence {e(t)} to be identically null. We assign a fixed value to the sparsity of the error φ(θ 3 ) (expressed in terms of the number φ(θ 3 ) 0 of nonzero components in φ(θ 3 )) and then solve problem [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF] 100 times on different independent simulations of input-output data of length N = 100 each. This procedure is repeated for different values of φ(θ 3 ) 0 reported in Table 1. In this table we display for each given value of φ(θ 3 ) 0 , the percentage of successes in attempting to compute the solution of ( 12) by solving [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF]. Since the data (x(t), y(t)), t = 1, . . . , N , is generated from [START_REF] Sun | Switched Linear Systems: Control and Design[END_REF] with white noise as input and uniformly distributed discrete mode, it holds with overwhelming probability that the columns of X are in general position. Consequently, it can be reasonably assumed that ν 4 (X) is as small as 4 i.e., ν 4 (X) is equal to the dimension of x(t), see [START_REF] Sun | Switched Linear Systems: Control and Design[END_REF]. According to Theorem 11, if φ(θ 3 ) 0 is roughly less than N -ν4(X)

2

= 48, then equivalence holds between problems ( 12) and ( 18) and the unique solution to both of them is θ 3 . This can be verified from the results of Table 1 where we see that [START_REF] Lauer | A product-of-errors framework for linear hybrid identification problem[END_REF] effectively solves [START_REF] Elhamifar | Sparse subspace clustering[END_REF] successfully with a score of 100% over 100 trials (on randomly generated data) once φ(θ 3 ) 0 falls under 48.

φ(θ 3 ) 0 58% 55% 53% 50% 48% 45% # succ. 46% 76% 94% 99% 100% 100%
Table 1 Equivalence between ℓ0 and ℓ1 minimizations versus the sparsity of φ(θ3). The ℓ0 norm of φ(θ3) is expressed as a fraction of the nonzero entries over the total length of the vector φ(θ3).

Now we propose, in the same conditions as the first experiment, to solve the sparse optimization problem [START_REF] Elhamifar | Sparse subspace clustering[END_REF] with the reweighted ℓ 1 optimization technique described in Algorithm 1 [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. The related results are presented in Table 2. It turns out that the reweighted ℓ 1 optimization approach significantly improves the basis pursuit technique. When φ(θ 3 ) 0 starts getting much larger than φ(θ 1 ) 0 and φ(θ 2 ) 0 , Algorithm 1 may not keep on converging towards θ 3 any longer. Instead, it is likely to converge towards θ 1 or θ 2 since φ(θ 1 ) 0 and φ(θ 2 ) 0 decrease as φ(θ 3 ) 0 increases.

φ(θ 3 ) 0 58% 55% 53% 50% 48% 45% # succ. 94% 100% 100% 100% 100% 100%
Table 2 Approximation of ℓ0 by reweighted ℓ1 minimization versus the sparsity of φ(θ3). The ℓ0 norm of φ(θ3) is expressed as a fraction of the nonzero entries over the total length of the vector φ(θ3).

Identification of the PVs

The second objective is to test the statistical robustness of the identification algorithm. For this purpose, we use 100 different independent realizations of the input, the discrete state and the output noise (SNR=30 dB) to generate 100 data sequences of length N = 600 each. The identification algorithm (Algorithm 2 indeed) is then run on each of these different 100 data sequences. The userspecified parameters of Algorithm 1 and Algorithm 2 are set to ε = 0.1, Tol = 0.001 and Thres = 0.05. At each run, the first 300 points are used to identify a model and the whole sequence of length 600 is used to validate the estimated model, i.e., to verify its ability to reconstruct the system output from the true input and an estimated discrete state. This is evaluated with the criterion [START_REF] Ljung | System Identification Toolbox User's Guide[END_REF] FIT

= 1 - ŷ -y 2 y -ȳ1 N 2 × 100% (31)
which measures the fitting error between the true output sequence y and the estimated model output sequence ŷ.

In this formula, ȳ stands for the mean of the true output sequence and 1 N is an N -dimensional vector with all entries equal to one. The reader can refer to Figure 2 for an insight into the form of the input-output signals. For better visualization purpose, the size of the observation window in that figure has been shortened.

The number of submodels is given a priori. We start by assuming that the number of submodels is fed into the identification algorithm. We present in Table 3 the average values of the estimated PVs together with their standard deviations over 100 independent runs of the algorithm. Along with those results are provided, for comparison purpose, the PVs' estimates the standard least squares would yield if the discrete mode sequence were fully known. By comparing the averaged estimates θi s i=1 , of the PVs displayed in Table 3 to the true values θ i s i=1 given in [START_REF] Vidal | Recursive identification of switched ARX systems[END_REF], we can see that the proposed algorithm has effectively recovered the true PVs with a relatively good precision despite the presence of noise. Moreover, by judging from the standard deviations σ( θi ), we are prompted to conclude that the algorithm performs well on statistically independent realizations of the input-output data. Of course, when the data is noise-free, the parameters are exactly recovered by the algorithm.

In Figure 3 is represented the distribution of the FIT over 100 of the identification algorithm on independent input-output data. This plot shows that most of the runs of the algorithm yield a FIT greater than 90%. In fact 98% of the runs produce a FIT measure larger than 87% (which means 100% if there were no noise in the data) on both identification and validation data. It can therefore be concluded that 98% of the runs yield the correct PVs. Least squares average estimates if the discrete state were known.

Table 3 Comparison of the proposed identification algorithm to standard least squares (if the discrete state were known) over 100 independent runs: ε = 0.1, Tol = 0.001, Thres = 0.05. # successes 98% 90% 96% 87% FIT 97.9% 93% 83.8% 63.7% Table 4 Performance of our method in determining the correct number of submodels, under different proportions of noise. Here, # successes is the number of times the algorithm yields the correct number of submodels (which is equal to 3) over 100 runs using a fixed threshold Thres. The value of FIT displayed here is the average fit over 100 independent runs.

The number of submodels is unknown. We turn now to the case when the number s of submodels is not known a priori. In such a context, we need to estimate, as in Algorithm 2, the number of submodels together with the parameter vectors. As we have seen in Subsection 3.5, the ability of the algorithm to provide the correct number (in fact the minimum number) of submodels from noisy data depends heavily on the threshold Thres6 . And selecting a threshold in a noisy environment is known to be challenging in general. For different levels of noise, the number of submodels is estimated over 100 independent realizations of data. The results displayed in Table 4 reveal a rather good tendency of the proposed method in recovering the true number of submodels.

Conclusion

The paper discusses a new optimization-based method for the identification of linear switched systems. By exploiting some ideas from the field of sparse signal recovery, we have first formulated this challenging iden-tification problem as a combinatorial ℓ 0 minimization problem. As such however, the problem is still computationally intractable so that practical implementation may involve relaxing it into an ℓ 1 norm based program, which is convex and therefore solvable with classical and well documented tools. Some sufficient conditions are derived for the convex relaxation strategy to exactly recover the solution of the initial ℓ 0 norm problem. In practice, different conclusive experiments on simulated data tend to show that the performance of the convex relaxation technique can be boosted far more beyond the established theoretical conditions of equivalence. Hence, the method for extracting the parameter vectors associated with the different subsystems, proves to be both computationally simple and satisfactorily efficient. As future work, we may envision to further study the properties of the method by (i) studying the convergence of the reweighted ℓ 1 minimization problem described in Algorithm 1 in the context of hybrid system identification, (ii) looking for tight bounds on the estimation error in highly noisy conditions. It would also be interesting to look for a way of efficiently extending the proposed method to switched systems in which the noise has a specific structure (for example, when the linear subsystems are described with ARMAX models).

  ) with λ t ∈ {1, 2, 3} and θ 1 = -0.40 0.25 -0.15 0.08 ⊤ , θ 2 = 1.55 -0.58 -2.10 0.96 ⊤ , θ 3 = 1 -0.24 -0.65 0.30 ⊤ .

  Average estimates over 100 independent runs of the identification algorithm: ε = 0.1, Tol = 0.001, Thres = 0.05.

Figure 2 .

 2 Figure 2. 100 identification and 100 validation data with an SNR of 30 dB. The algorithm is run with ε = 0.1, Tol = 0.001 and Thres = 0.05: Obtained FIT = 93%.

Figure 3 .

 3 Figure 3. Distribution of the FIT measure over the 100 independent runs of the algorithm. ε = 0.1, Tol = 0.001, Thres = 0.05.

Strictly speaking, ℓ0 is not a norm as it does not satisfy the property of positive scalability, i.e., λz 0 = |λ| z 0 does not hold in general.

We say that a vector x is sparser than another vector y of the same dimension if x 0 ≤ y 0 .

This problem is common to most existing methods that estimate the number of submodels, i.e. their capacity of determining the correct number of submodels depends on a user-specified parameter, see[START_REF] Bemporad | A bounded-error approach to piecewise affine system identification[END_REF][START_REF] Bako | Algebraic identification of switched MIMO ARX models[END_REF][START_REF] Boukharouba | Identification of piecewise affine systems based on dempster-shafer theory[END_REF][START_REF] Ozay | A sparsification approach to set membership identification of a class of affine hybrid systems[END_REF].

Acknowledgements

The author would like to thank the anonymous reviewers, the editor and the associate editor whose constructive comments on an earlier version of this manuscript have been of great help in improving the presentation.


Kiran Kumar Katakam

Peter Nejsum

Niels Christian Kyvsgaard

Claus B Jørgensen

Stig Milan Thamsborg

Molecular and parasitological tools for the study of Ascaridia galli population dynamics in chickens

Keywords: Ascaridia galli, larvae recovery method, cohorts, population dynamics

Molecular and parasitological tools for the study of Ascaridia galli population dynamics in chickens

Abstract

Experiments were first conducted to compare and evaluate different methods of Ascaridia galli larval recovery from the chicken intestine. The number of larvae recovered from the intestinal wall of chickens infected with 1000 embryonated A. galli eggs and killed 15 days post infection (p.i.) by three methods (EDTA, pepsin digestion and scraping) were compared. The EDTA and pepsin digestion were found to be the most efficient methods with no significant difference (P > 0.05) in number of recovered larvae between the two. Subsequently, three different A. galli cohorts were established using the polymerase chain reaction linked restriction fragment length polymorphism (PCR-RFLP) technique. A 533 bp long region of the cytochrome c oxidase subunit 1 (cox1) gene of the mitochondrial DNA (mtDNA) was targeted and 22 A. galli females were allocated in to three different haplotypes. The four females with the highest embryonation rate from each halpotype group (total 12 females) were selected and used to inoculate each of 12 chickens with a dose of 1000 embryonated eggs. The chickens were killed 15 days p.i. and A. galli larvae were recovered from the small intestinal wall by the EDTA method and by sieving the lumen content on a 90 µm sieve. DNA of 40 larvae from each of the three different haplotypes was extracted using a worm lysis buffer (WLB) and PCR-RFLP analysis of these larvae revealed same haplotype as that of their maternal parent. The identification of distinguishable cohorts may be a powerful tool in population studies of parasite turnover within the animal host.

Introduction

Ascaridia galli (Schrank, 1788) is a nematode parasite that infects domestic as well as wild birds [START_REF] Soulsby | Helminths, Arthropods and Protozoa of Domesticated Animals, 7th edition[END_REF]. Studies have suggested that A. galli is the most common nematode in all types of production systems and has a worldwide distribution. [START_REF] Permin | On investigation in different media for embryonation of Ascaridia galli eggs[END_REF][START_REF] Ashenafi | Study on gastrointestinal helminths of local chickens in central Ethiopia[END_REF]Martin-Pacho et al., 2005;[START_REF] Rabbi | Gastrointestinal helminths infection in different types of poultry. Bangl[END_REF][START_REF] Abdelqader | Prevalence and burden of gastrointestinal helminthes among local chickens, in northern Jordan[END_REF]. The adult stage inhabits the lumen of the small intestine of the host, feeding mostly on ingesta [START_REF] Ackert | The food of the fowl nematode, Ascaridia lineata (Schneider)[END_REF][START_REF] Soulsby | Helminths, Arthropods and Protozoa of Domesticated Animals, 7th edition[END_REF]. A. galli infections have been associated with reductions in egg production in laying hens and in overall growth in chickens [START_REF] Soulsby | Helminths, Arthropods and Protozoa of Domesticated Animals, 7th edition[END_REF][START_REF] Ramadan | Some pathological and biochemical studies on experimental ascaridiasis in chickens[END_REF]. A. galli may also play a role in transmission of Salmonella infections [START_REF] Chadfield | Investigation of the parasitic nematode Ascaridia galli (Shrank 1788) as a potential vector for Salmonella eneterica dissemination in poultry[END_REF][START_REF] Eigaard | Infection and excretion of Salmonella enteritidis in two different chicken lines with concurrent Ascaridia galli infection[END_REF]) and avian reoviruses [START_REF] Calnek | Diseases of Poultry[END_REF] resulting in disease and economic loses.

The life cycle of A. galli is direct and the infective stage is eggs harbouring third stage larvae. After ingestion, eggs reach the duodenum and hatch within 24 hours [START_REF] Ackert | On the habitat of Ascaridia perspicillum (Rud)[END_REF], and larvae are released into the lumen of the intestine. The larvae then enter the mucosa of the small intestine and most larvae initiate a histotrophic phase between day 8 to 17 post infection (p.i.) [START_REF] Tugwell | On the tissue phase of the life cycle of the fowl nematode Ascaridia galli (Schrank)[END_REF]Ackert 1952, Herd and[START_REF] Herd | Arrested development and the histotropic phase of Ascaridia galli in the chicken[END_REF]. Both onset and length of this phase depend on infection dose [START_REF] Herd | Arrested development and the histotropic phase of Ascaridia galli in the chicken[END_REF].

New regulations in the European Union, setting standards for protection and welfare of laying hens (Anon, 1999), are implemented from 1st January 2012 and prescribe substitution of traditional cage systems with floor systems and enriched cage systems [START_REF] Willer | The world of organic agriculture: Statistics and emerging trends, International Federation of Organic Agriculture Movements (IFOAM)[END_REF].

In floor and free-range systems, the risk of A. galli infection is known to be very high [START_REF] Permin | On investigation in different media for embryonation of Ascaridia galli eggs[END_REF]. In addition, the increased demand for more natural or ethical responsible products in recent years is making organic farming systems more common. As a consequence, the overall Understanding the host-parasite relationships is a necessary prerequisite to develop such control strategies. For this purpose tools for recovery of the larval stages and identification of subpopulations will be very useful. The present study aimed at developing appropriate research tools for understanding the population dynamics of A. galli. A standard larval recovery technique for post mortem examination and a method to establish genetically different A. galli cohorts were therefore developed and evaluated.

Materials and Methods

Preparation of egg batches. 30 Hy-line breed chickens collected at an organic farm (no use of anthelmintics for 10 years) were killed by cervical disarticulation, eviscerated and a total of 30 adult A. galli female worms were collected from 20 hosts. A. galli eggs were isolated from the uterus of these worms separately and incubated for 6 weeks in 0.1N sulphuric acid [START_REF] Permin | On investigation in different media for embryonation of Ascaridia galli eggs[END_REF].

Twenty two egg batches were selected based on the highest embryonation rates (above 50%).

Experiment 1: Larval recovery. Infection of chickens. Twenty eight four-week old Lohmann Silver chickens were obtained from a conventional indoor pullet raising farm, supposedly uninfected premises, and randomly allocated to three groups of 10, 9 and 9 chickens and housed in separate sheds. After 2 days of acclimatization, 7 chickens in each shed were orally infected with a dose of 1000 embryonated eggs and the rest kept as uninfected controls. The egg batch was made as a mixture of eggs originating from the 22 females mentioned above in equal proportions. Processing the intestinal contents. 15 days p.i. 27 chickens including 6 controls were killed by cervical disarticulation and randomly allocated to 3 larval recovery methods (see below) after stratification for shed. Pepsin digestion method was applied for controls. The gastrointestinal tract was removed from gizzard to cloaca. The intestines were cut open longitudinally and washed gently. All the luminal contents including the washings were collected on a 90 µm sieve and washed thoroughly with tap water. The retained material was transferred to a 50 ml tube and stored in 70 % ethanol at 5°C. The intestinal wall was subsequently subjected to either of three different methods: 1) EDTA, 2) pepsin digestion, 3) scraping. Furthermore, after initial processing by EDTA or scraping, intestines were pepsin digested.

EDTA method. 10mM EDTA in 0.9% saline solution was prepared [START_REF] Kringel | A method for the recovery of immature Trichuris suis from pig intestine[END_REF]. The small intestine was cut into 4 pieces longitudinally and hung in 10 mM EDTA in 0.9% NaCl for 6 hours in a 500 ml conical glass. After removing the intestinal wall the contents were allowed to settle for 30 minutes. The supernatant was discarded and the sediment was collected and stored in 70% ethanol at 5°C. Pepsin digestion method. The digestion fluid of 1L was prepared freshly [START_REF] Kapel | Infectivity, persistence, and antibody response to domestic and sylvatic Trichinella spp. In experimentally infected pigs[END_REF].

The small intestine was cut in to 2-3 cm pieces and digested in a beaker with 100 ml digestion fluid under moderate stirring (125 rpm) at 38°C for 45 min. The undigested intestinal pieces were rinsed on a metal sieve of 2-3 mm size with water and digested material was transferred to a conical glass, filled with water and left for 20 min. The supernatant was poured off and again filled with water and allowed to sediment for 20 minutes. This procedure was repeated for 5-6 times. The sediment was collected and stored in 70% ethanol at 5°C. Scraping method. Using a glass slide the mucosal and sub-mucosal layers were carefully scraped off the small intestine and collected on a 90 µm sieve and washed thoroughly with water. The retained material was collected and stored in 70% ethanol at 5°C.

Measurement of length of A. galli larvae.

With the help of tweezers, the length of each larva was measured manually using a measuring scale.

Comparison of quality of DNA of larvae following recovery. DNA from two times 20 larvae recovered by the EDTA and pepsin digestion methods, respectively, was extracted as described in the following section, PCR was conducted and bands were visualised after gel electrophoresis.

Experiment 2: Cohort study. Extraction of DNA. DNA was extracted from the 22 selected adult worms using the MasterPure™ DNA purification kit (Epicentre Biotechnologies) according to the manufacturer's protocol. DNA from single larvae was extracted using a worm lysis buffer [START_REF] Maafi | Molecular identification of cyst-forming nematodes (Heteroderidae) from Iran and a phylogeny based on ITS-rDNA sequences[END_REF], and carried out according to [START_REF] Nejsum | A novel technique for identification of Ascaris suum cohorts in pigs[END_REF].

Selection of the target region (cox1 gene) in mtDNA. Sequence information of the cox1 gene of the mtDNA was obtained from Ascaris suum, Caenorhabditis elegans and Necator americanus using BLAST (http://www.ncbi.nlm.nih.gov/blast/Blast.cgi). Conserved regions were subsequently identified by aligning the sequences using TCoffee (http://igsserver.cnrsmrs.fr/Tcoffee/tcoffee_cgi/index.cgi) and used as target regions for primers. The forward (GCox18F-5': TTGTGGTCTGGTATGGTTGG) and reverse (GCox18R-5: PCR and sequencing. The PCR was performed in a TECHNE ® TC-512 thermal cycler. The thermal profile for the PCR reaction was: 15 minutes at 95°C for initial activation of Taq polymerase followed by 35 cycles consisting of: 30 seconds at 95°C (denaturation), 40 seconds at 55°C (annealing), 1 minute at 72°C (extension) and a final elongation of 10 minutes at 72°C.

Negative H 2 O controls were included in all runs. The amplified fragment was purified by Qiaquick PCR-purification kit (Qiagen, Germany) and sequenced in both directions using a BigDye terminator sequencing system according to the manufacturer's instructions (Applied Biosystems, USA). After ethanol precipitation the sequencing products were run on an ABI3130XL (Applied Biosystems, USA).

Design of new primers. The sequence information was used to design new primers (GCox1F4F: 5'-ATTATTACTGCTCATGCTATTTTGATG and GCox14R: 5'-CAAAACAAATGTTGATAAATCAAA GG) in conserved regions of the A. galli sequence giving an amplicon with a length of 533 bp.

Selection of restriction enzymes. After obtaining the sequence information, single nucleotide polymorphisms (SNP) were identified manually and NEB cutter V2.0 (http://tools.neb.com/NEBcutter2/index.php) was used to identify two restriction enzymes (NlaIII and SfcI) that could allocate the 22 samples into three haplotypic groups and also had the property of sharing similar buffers and reaction conditions. The PCR products of all 22 A. galli samples were digested overnight at 37°C using 3 units of both NlaIII and SfcI (New England, Biolabs, USA) and Infection of chicken. Twenty-five nine-week old chickens of Lohmann Sandy hybrid were obtained from another farm than in the former experiment. Ten chickens were immediately killed and all found to be A. galli negative. The remaining 15 chickens were randomly allocated into three groups of 5 and housed in three different sheds. Four chickens in each group were inoculated with 1000 embryonated A. galli eggs of the three different haplotypes and one chicken in each shed was kept as uninfected control. 15 days p.i. all chickens were killed and larvae recovered from the intestines by the EDTA method.

PCR-RFLP on larvae. The DNA of both the selected adult worms (positive controls) and larvae was subjected to PCR-RFLP as described above.

Data analysis. The recoveries of larvae from intestinal lumen and wall by different methods were compared by a generalized linear model (Proc Genmod, SAS v 9.1) by specifying a negative binomial distribution of larval counts. The difference in the length of the larvae recovered by different methods from the intestinal wall was analysed in a mixed effect variance model with the method as the fixed effect and the chicken as a random effect. The difference in the length of larvae between intestinal content and intestinal wall was analysed in a mixed effect variance model with location (lumen or wall) as the fixed effect and the chicken as the random effect.

One-way analysis of variance (ANOVA) was used to compare the recovery of larvae in the intestinal lumen, wall and total respectively, of three different haplotypes in the cohort experiment.

Animal ethics

The study was carried out in accordance with the requirements of The Danish Animal Ethics Committee (permit # 2005/561-1060)

Results

Experiment 1: Larval recovery. The percentage of total larvae recovered by different methods ranged from 0.83-2.56% of the infection dose (Table 1). The recoveries from intestinal wall were overall significantly different (P = 0.018) and subsequent pair wise analysis showed significant higher recovery by the EDTA method compared with scraping (P = 0.006) and higher recovery by pepsin digestion than scraping (P = 0.003). Pepsin digestion of the intestines after EDTA incubation yielded no more larvae while digestion after scraping yielded a single larva.

The recovered larvae had an overall mean length of 3.1mm. There was no statistical difference in the mean length of larvae recovered from the intestinal wall by either EDTA method or pepsin digestion method (P = 0.66). Similarly there was no significant difference in mean length of larvae recovered from the lumen and those recovered from the intestinal wall (P = 0.90).

Three chickens harboured nine adult worms indicating previous exposure to A. galli.

For all larvae (n = 40) recovered by EDTA and pepsin digestion, it was possible to amplify the cox1 gene by PCR.

Experiment 2: Cohort study. In the three groups infected with eggs of haplotype I, II and III respectively, a mean percentage of 2.2, 2.6 and 1.2 of the infection dose with a mean number of (individual numbers 8, 4, 21 and 15) larvae were recovered. There was no significant difference in larval recovery between different haplotype groups. Forty larvae randomly selected from all 4 chickens of each of three maternal-haplotype groups (total 120) recovered from the chickens were PCR-RFLP typed and all were found to hold the same band pattern as their female of origin. The band pattern of some of the larvae is shown in Figure 1. The nucleotide sequences for the three haplotypes have been submitted to the GenBank under the following accession numbers GU138668-GU138670.

Discussion

Previous studies reported that the histotrophic phase is a normal part of A. galli life cycle [START_REF] Herd | Arrested development and the histotropic phase of Ascaridia galli in the chicken[END_REF] and that most histotrophic larvae are in the mucosal layer and only rarely penetrate deeper into the intestinal wall [START_REF] Ackert | On the habitat of Ascaridia perspicillum (Rud)[END_REF]. A commonly used method for isolation of histotrophic larvae has been pepsin digestion e.g. as used by Herd and Mc Naught (1975) who incubated intestinal wall for 4 -12 h at 37°C. However, [START_REF] Irvine | Efficient polymerase chain reaction detection of the second internal transcribed spacer of mucosa-derived larvae is dependent on the larval extraction method[END_REF] reported that digestion for prolonged period (more than 2h) may be associated with reduced PCR efficiency. We therefore compared a pepsin digestion technique developed for the release of Trichinella spp. larvae [START_REF] Kapel | Infectivity, persistence, and antibody response to domestic and sylvatic Trichinella spp. In experimentally infected pigs[END_REF] but with reduced digestion period (45 minutes) with the EDTA method which has been developed for the recovery of small immature Trichuris suis larvae from large intestine of pigs [START_REF] Kringel | A method for the recovery of immature Trichuris suis from pig intestine[END_REF]. These methods were compared with manual scraping. Our study showed no difference in larval recovery by the EDTA and pepsin digestion methods and that both were superior compared to scraping. Following scraping it was very difficult to find the larvae embedded in the clumps of mucosa. EDTA presumably loosens the connections between mucosal cells whereby embedded larvae are released whereas digestion with acid pepsin dissolves the mucosa and releases the larvae. The fact that A. galli larvae rarely penetrate deep in to intestinal tissue [START_REF] Ackert | On the habitat of Ascaridia perspicillum (Rud)[END_REF]) might be the reason for not finding any additional larvae when the intestines were pepsin digested after EDTA incubation. Based on this, and the fact that these samples were the easiest samples to read, the EDTA method is recommended for quantitative recovery.

We found no difference with respect to length neither between lumen and wall nor between EDTA and pepsin digestion method. It might have been suspected that a given larval recovery method would be better to recover larvae of a certain size, but we did not find any indication for such a selection in our study although the larvae had a relatively large variation in size.

Applying PCR-RFLP on the mitochondrial DNA it was possible to establish three different types of cohorts. We haplotyped 120 larvae and all showed a similar band pattern on the gel as that of their maternal parent. This supports the underlying hypothesis that mtDNA is maternally inherited in A. galli as is the case for e.g. Ascaris suum [START_REF] Nejsum | A novel technique for identification of Ascaris suum cohorts in pigs[END_REF]) and mammals in general [START_REF] Hutchison | Maternal inheritance of mammalian mitochondrial DNA[END_REF].

The possibility to identify different cohorts makes it possible to follow different groups of A. galli, their location in the host and their relative establishment rates. [START_REF] Dobson | Population-dynamics of Trichostrongylus colubriformis in sheep-the effect of infection rate on the establishment of infective larvae and parasite fecundity[END_REF] made a cohort study using levamisole resistant and susceptible strains of Trichostrongylus colubriformis whereas [START_REF] Roepstorff | Experimental Oesophagostomum dentatum infections in the pig: worm populations resulting from trickle infections with three dose levels of larvae[END_REF] used pyrantel resistant and susceptible strains of Oesophagostomum dentatum in his study. But this approach to cohort studies results only in two different cohorts and depends on availability of solidly resistant strains. The latter is not the case for a range of nematodes. In contrast [START_REF] Sørensen | Elucidation of Schistosoma japonicum population dynamics in pigs using PCR-based identification of individuals representing distinct cohorts[END_REF] PCR-RFLP method on the mtDNA and were able to establish four different cohorts of A. suum. The utilization of the natural genetic variation in the mtDNA to identify different cohorts of parasites therefore not only makes it possible to identify several different cohorts but also an ability to work with naturally occurring worms in contrast to isolates or specific strains. However, similar establishment rates need to be confirmed before embarking on comparisons.

The development of reliable quantitative recovery methods in combination with genetically marked cohorts provides powerful tools for the study of population dynamics of nematodes in poultry and other livestock.

F o r P e e r R e v i e w O n l y

  of A. galli is expected to increase in coming years, and alternative control strategies are warranted in order to avoid increasing dependence on repeated use of anthelmentics.

  which were designed using Primer3 (http://frodo.wi.mit.edu/) and flanked a region of ~800 bp.

  were subjected to gel electrophoresis on a 1.5% agarose gel (100 V, 90 minutes).

  individual numbers 1, 28, 43 and 14), 25.7 (individual numbers 3, 7, 85 and 8) and 12.0

  used PCR-RFLP technique to distinguish between two isolates of Schistosoma japonicum whereas[START_REF] Nejsum | A novel technique for identification of Ascaris suum cohorts in pigs[END_REF] used the present

Figure 1 .

 1 Figure 1. PCR-RFLP on mtDNA. M: 100 bp size marker, lane 1 and 2-4: band patterns of adult

Table 1 .

 1 The number of A. galli larvae recovered from intestinal wall and lumen contents of chicken experimentally infected with 1000 embryonated eggs and euthanized 15 days p.i.Intestinal wall recovery by EDTA, pepsin digestion or scraping methods.

		Intestinal wall		Lumen contents
		Mean length		Mean length
	Number of larvae	(mm)	Number of larvae	(mm)	Total larvae
	EDTA method
	Mean (±s.d.)	5.4(±8.2)	2.5(±3.2)	21.3(±20.8)	4.4(±3)	26.7(±23.6)
	Range	0-23	1-9.5	3-54	1.5-7	5-62
	Mean recovery (%)	0.5%		2.1%		2.7%
	Pepsin digestion method
	Mean (±s.d.)	6.7(±10.5)	2.5±1.6	20.1(±16.4)	3±1.4	26.7 (±20.9)
	Range	0-30	1.5-4	2-47	1-8.5	3-57
	Mean recovery (%)	0.7 %		2.0%		2.7 %
	Scraping method
	Mean(±s.d.)	0.6(±0.53)	2.6±0.2	7.6(±9.84)	3.2±1.2	8.9(±10.3)
	Range	0-1	2.5-3	0-25	1.5-10	0-26
	Mean recovery(%)	0.06%		0.8%		0.9%

172x92mm (96 x 96 DPI)

E-mail: cavanagh@metronet.co.uk URL: http://mc.manuscriptcentral.com/cavp Avian Pathology

Acknowledgements

The Danish Centre for Experimental Parasitology is acknowledged for this project. Thanks to Minna Jakobsen for technical assistance on DNA sequencing and Lise-Lotte Christiansen for assistance with the larval recovery experiment.