Thibault Defourneau

Mihaly Petreczky

Realization theory of recurrent neural networks and rational systems

In this paper, we show that, under mild assumptions, input-output behavior of a continuous-time recurrent neural network (RNN) can be represented by a rational or polynomial non-linear system. The assumptions concern the activation function of RNNs, and they are satisfied by many classical activation functions. We also present an algorithm for constructing the polynomial and rational system. This embedding of RNNs into rational systems can be useful for stability, identifiability, and realization theory for RNNs, as these problems have been studied for polynomial/rational systems. In particular, we use this embedding for deriving necessary conditions for realizability of an input-output map by RNN, and for deriving sufficient conditions for minimality of an RNN.

I. INTRODUCTION

One of the challenges in machine learning is to provide a mathematical theory for analyzing learning algorithms. Recently, there has been a surge of interest in the use of neural networks, leading to the emergence of the field of deep learning. One of the most widespread models used in deep learning are recurrent neural networks (RNNs). RNNs can seen as non-linear dynamical systems equipped with an internal state, input and output. Learning such an RNN from data is equivalent to estimating the parameters of the RNN, viewed as a dynamical system. That is, learning algorithms for RNNs correspond to system identification algorithms, and developing a mathematical theory for learning RNNs is equivalent to developing system identification for RNNs. There is a rich literature on system identification, in particular on system identification for linear systems [START_REF] Ljung | System identification: theory for the user[END_REF]. Note that linear dynamical systems are a particular class of RNNs.

One of the principal building blocks of system identification theory for linear systems is realization theory. Realization theory can be viewed as an attempt to solve an idealized system identification problem, where there is infinite data, not modelling error, etc. In general, the aim of realization theory is to understand the relationship between an observed behavior and dynamical systems producing this observed behavior. For the particular case of RNNs, the main questions of realization theory can be stated as follows:

1) Which class of observed behaviors (input-output maps) can be represented by an RNN ?

Centre de Recherche en Informatique, Signal et Automatique de Lille (CRIStAL), CNRS, Ecole Centrale, Universit de Lille, Avenue Carl Gauss, 59650 Villeneuve-d'Ascq, France (e-mail:thibault.defourneau@univ-lille.fr, mihaly.petreczky@centralelille.fr)

This work was partially funded by: (1) CPER Data project, which is co-financed by European Union with the financial support of European Regional Development Fund (ERDF), French State and the French Region of Hauts-de-Francen, (2) by CNRS project PEPS blanc BayesRealForRNN, and (3) by Agence Nationale pour la Recherche Project ROCC-SYS under the Grant agreement ANR-14-C27-0008

2) How can we characterize minimal RNNs (RNNs of the least complexity) representing a certain observed behavior ? What is the appropriate definition of minimality (smallest number of neurons, etc.) for RNNs, are minimal RNNs are related by some transformation ? 3) Is there a constructive procedure for constructing a RNN representation from input-output behavior which can be proven to be mathematically correct ? For linear systems, realization theory [START_REF] Kailath | Linear Systems[END_REF], [START_REF] Lindquist | Linear Stochastic Systems[END_REF] has been useful for system identification, for example it helped to address identifiability, canonical forms and gave rise to subspace identification algorithms. We expect that realization theory of RNNs will lead to similarly useful results for the latter dynamical systems.

In order to develop realization theory of RNNs, we embed RNNs into the class of rational systems, and then we use realization theory of rational/polynomial systems in order to derive new results on realization theory of RNNs.

In this paper, we consider RNNs in continuous-time, in order to avoid some technical difficulties encountered in the discrete-time case. Note that discrete-time RNNs can be viewed as discretizations of continuous-time RNNs.

In this paper we assume that the high-order derivatives of the activation function satisfy a polynomial equation. Several widely-used activation functions have this property.

• We show that an input-output map can be realized by a RNN, only if it can be realized by a rational system, i.e. a non-linear system defined by vector fields and readout maps which are fractions of polynomials. We present an explicit construction of such a rational system. • We present a necessary condition for existence of a realization by RNNs, using results from realization theory of rational systems. This necessary condition is a generalization of the well-known rank condition for Hankel matrices of linear systems. • We formulate sufficient conditions for observability/reachability/minimality of RNNs, using existing realization theory for rational systems [START_REF] Němcová | Realization theory for rational systems: The existence of rational realizations[END_REF], [START_REF] Němcová | Rational Systems in Control and System Theory[END_REF], [START_REF] Bartoszewicz | Minimal polynomial realizations[END_REF], [START_REF] Wang | Algebraic differential equations and rational control systems[END_REF]. Note that RNNs could be viewed as analytic systems and one could try to apply realization theory of analytic systems [START_REF] Jakubczyk | Realization theory for nonlinear systems: three approaches[END_REF], [START_REF] Isidori | Nonlinear control systems[END_REF], [START_REF] Hermann | Nonlinear controllability and observability[END_REF]. However, realization theory for analytic systems is not computationally effective. In contrast, the conditions for observability/reachability/minimality of RNNs which are derived in this paper are computational effective and they are less restrictive than those which can be obtained by viewing RNNs as analytic systems [START_REF] Jakubczyk | Realization theory for nonlinear systems: three approaches[END_REF], [START_REF] Isidori | Nonlinear control systems[END_REF], [START_REF] Hermann | Nonlinear controllability and observability[END_REF].

To the best of our knowledge, the results of the paper are new. RNNs have been widely used in the machine learning literature, both in discrete-time and continuoustime, [START_REF] Sustkever | Training Recurrent Neural Networks[END_REF], [START_REF] Caterini | Deep Neural Networks in a Mathematical Framework[END_REF]. Observability of RNNs was studied in [START_REF] Albertini | State observability in recurrent neural networks[END_REF], controllability in [START_REF] Qiao | Further results on controllability of recurrent neural networks[END_REF] and minimality in [START_REF] Albertini | For Neural Networks, Function Determines Form[END_REF]. This paper was inspired by [START_REF] Albertini | State observability in recurrent neural networks[END_REF], [START_REF] Qiao | Further results on controllability of recurrent neural networks[END_REF], [START_REF] Albertini | For Neural Networks, Function Determines Form[END_REF], but in contrast to [START_REF] Albertini | State observability in recurrent neural networks[END_REF], [START_REF] Qiao | Further results on controllability of recurrent neural networks[END_REF], [START_REF] Albertini | For Neural Networks, Function Determines Form[END_REF], we do not use any assumption on the structure of the weight matrices, except for observability issues. This means that the results of this paper can be applied even when the results of [START_REF] Albertini | State observability in recurrent neural networks[END_REF], [START_REF] Qiao | Further results on controllability of recurrent neural networks[END_REF], [START_REF] Albertini | For Neural Networks, Function Determines Form[END_REF] are not applicable.

In Section II we present the basic notation and terminology. In Section III we present the construction of the rational system which realizes the same input-output map as an RNN. In Section IV, we use the results of Section III to derive necessary conditions for existence of a realization by RNN. Finally, in Section V we use the results of Section III to derive sufficient conditions for minimality of RNNs. An extended version of this paper containing detailed proofs can be found in the technical report [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

II. BASIC DEFINITIONS

In this section, we fix the notation used in the paper and we recall some algebraic tools necessary for this paper. Then we recall the definition of RNNs and of rational systems.

A. Preliminaries

We use the standard terminology and notation from commutative algebra and algebraic geometry see [START_REF] Kunz | Introduction to commutative algebra and algebraic geometry[END_REF], [START_REF] Zariski | Commutative algebra I[END_REF], [START_REF] Cox | Ideals, varieties, and algorithms: An introduction to computational algebraic geometry and commutative algebra[END_REF]. In particular, by R[X 1 , . . . , X n ] we denote the algebra of real polynomials in n variables and by R(X 1 , . . . , X n ) we denote the quotient field of R[X 1 , . . . , X n ]. We refer to the element of R(X 1 , . . . , X n ) as rational functions in n variables.

B. Recurrent neural networks.

Below we define formally what we mean by recurrent neural networks in continuous-time. We will follow the notation of [START_REF] Albertini | State observability in recurrent neural networks[END_REF], [START_REF] Qiao | Further results on controllability of recurrent neural networks[END_REF], [START_REF] Albertini | For Neural Networks, Function Determines Form[END_REF].

Definition 1: A recurrent neural network, abbreviated by RNN, with input-space U ⊂ R m and output-space R p , is a dynamical system

Σ : ẋ(t) = - → σ Ax(t) + Bu(t) , x(0) = x 0 y(t) = Cx(t) (1)
where

• σ : R → R is a continuous globally Lipschitz scalar function. It is called the activation function. • A ∈ R n×n , B ∈ R n×m and C ∈ R p×n are matrices, called weight matrices, • the map - → σ : R n → R n is defined by - → σ : (x 1 , . . . , x n ) T → (σ(x 1 ), . . . , σ(x n )) T ,
• u(t) is an input, x(t) ∈ R n is the state and y(t) is the output at time t, and x 0 ∈ R n is the initial-state. We denote such a system by Σ = (A, B, C, U, σ, x 0 ).

If we take σ the identity map in Definition 1, it is the same as a linear system in control theory. Next we define formally what we mean by a solution and an input-output map of a system Σ = (A, B, C, U, σ, x 0 ). To this end, in the sequel, we denote by

P C([0; +∞[, X) the set of piecewise- continuous functions from [0; +∞[ to X ⊆ R k , k > 0. Definition 2: A triple (x, u, y) is a solution of an RNN Σ = (A, B, C, U, σ, x 0 ) if u ∈ P C([0; +∞[, U), x : [0; +∞[→ R n , y : [0; +∞[→ R p ,
x is absolutely continuous, and (1) holds.

Remark 1: As the activation function σ is globally Lipschitz, it follows that for every u ∈ P C([0; +∞[, U), there exist unique x, y such that (x, u, y) is a solution of Σ.

In this paper, we focus on solutions (x, u, y) of an RNN such that u is piecewise constant.

Notation 1 (U pc ): We denote by U pc the set of all piecewise-constant functions of the form u : [0; +∞[→ U.

Definition 3: Let p : U pc → P C([0; +∞[, R p ) be an input-output map, and let Σ = (A, B, C, U, σ, x 0 ) be an RNN. Σ is said to be a realization of the input-output map p if for every u ∈ U pc , the unique solution (x, u, y), x(0) = x 0 of Σ is such that p(u) = y.

C. Rational and polynomial systems.

Informally, rational respectively polynomial systems are control systems in continuous time, whose differential equations and readout maps are rational functions, i.e. they are fractions of polynomials, respectively polynomial functions.

Definition 4 (Polynomial and rational systems): A rational system with input-space U ⊂ R m , state-space R n and output-space R p is a dynamical system as

R : υi (t) = P i,u(t) (υ(t)) Q i,u(t) (υ(t)) , i = 1, . . . , n , υ(0) = υ 0 y k (t) = h k,1 (υ(t)) h k,2 (υ(t))) , k = 1, . . . , p (2)
where • u(t) is the input, υ(t) = (υ 1 (t), . . . , υ n (t)) T is the state and y(t) = (y 1 (t), . . . , y p (t)) T is the output at time t.

Moreover υ 0 ∈ R n is the initial state; • h k,1 , h k,2 , k = 1, . . . , p are non-zero polynomials in n
variables and for all α ∈ U, i = 1, . . . , n, P i,α , Q i,α are polynomials in n variables, Q i,α is non-zero. We will identify the rational system R with the tuple

({P i,α , Q i,α } i=1,...,n,α∈U , {h k,1 , h k,2 } p k=1 , U, υ 0 ). We will say that R is polynomial, if h k,2 = 1, Q i,α = 1 for all k = 1, . . . , p, i = 1, . . . , n, α ∈ U.
Note that the existence and uniqueness of a state trajectory of a rational system requires some care, as in our definition we did not exclude the possibility that x(t) passes through the zero set of a denominator Q i,α .

Definition 5: A triplet (υ, u, y) is a solution of a rational system R of the form (2), if the input u : [0; +∞[→ U is piecewise constant, the state υ : [0; +∞[→ R n is absolutely continuous, the output y : [0; +∞[→ R p is piecewise continuous, and they satisfy

υi (t) Q i,u(t) (υ(t)) = P i,u(t) (υ(t)) y k (t) h k,2 (υ(t)) = h k,1 (υ(t)) , (3)
for 1 i n, and 1 k p, where υ(t) = (υ 1 (t), . . . , υ n (t)) T is the state and y(t) = (y 1 (t), . . . , y p (t)) T .

By uniqueness of a solution of an analytic differential equation, for any initial state υ 0 ∈ R n such that Q i,u(t) (υ 0 ) = 0 and h k,1 (υ 0 ) = 0, there exist at most one solution (υ, u, y) of R such that υ(0) = υ 0 . In particular, if R is a polynomial system, then for any initial state υ 0 ∈ R n there exist at most one solution (υ, u, y) of R such that υ(0) = υ 0 .

Definition 6: Let p : U pc → P C([0; +∞[, R p ) be an input-output map, and let R be a rational system of the form [START_REF] Albertini | For Neural Networks, Function Determines Form[END_REF]. We say that R realizes p, if for any u ∈ U pc there exists a solution (υ, u, y) of R such that υ(0) = υ 0 and p(u) = y.

III. REPRESENTING RNNS BY RATIONAL SYSTEMS.

In this section we show that an RNN realization of a given input-output map imply the existence of a rational system which is a realization of the same map. Moreover, we present the construction of such a rational realization. We introduce the following assumption.

Assumption 1 (A1): The function σ : R → R is analytic, and there exist an integer N > 0 and N analytic functions ξ 1 , . . . , ξ N : R → R such that

   σ V 0 (ξ 1 , . . . , ξ N ) = U 0 (ξ 1 , . . . , ξ N ) ξi V i (ξ 1 , . . . , ξ N ) = U i (ξ 1 , . . . , ξ N ) , if 1 i N (4)
where U k , V k are polynomials in N variables.

Assumption (A1) involves existence of analytic functions {ξ i } N i=1 and hence it is not easy to check it. In fact, Assumption (A1) can be replaced by the following hypothesis, which involves only derivatives of the activation function σ.

Assumption 2 (A2): The function σ : R → R is analytic, and there exist an integer N > 0 and a no-zero polynomial Q in N + 1 variables, such that Q(σ, σ (1) . . . , σ (N ) ) = 0 [START_REF] Bartoszewicz | Rational systems and observation fields[END_REF] where σ (i) denotes the i-th derivative of σ. Lemma 1 (Equivalence of (A1) and (A2)): A function σ satisfies Assumption (A1) if and only if it satisfies Assumption (A2). The proof of Lemma 1 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

Example 1: Let consider an analytic RNN of a response map for which the activation function σ is the hyperbolic tangent th, or the sigmoid function S given below:

∀x ∈ R , th(x) = e x -e -x e x + e -x , S(x) = 1 1 + e -x .
These functions are analytic and satisfy a differential polynomial equation, more precisely the hyperbolic tangent verifies y (1) = 1 -y 2 with y(0) = 0, and the sigmoid function satisfies y (1) = y(1 -y) with y(0) = 1 2 . It follows that Assumption (A2) holds.

Then we restrict the set of input maps, by supposing the following assumption.

Assumption 3 (Finite input set): In the rest of the paper we assume that U ⊂ R m is a finite set.

Notation 2: We denote by |U| the cardinality of U. We set U = {α 1 , . . . , α K }, where α i ∈ R m , and α i = α j if i = j. In that case, we have |U| = K.

The assumption that U is finite is not an overly restrictive one, and it is satisfied in many potential applications.

Next we present the main theorem of this paper. Theorem 1 (Embedding RNNs into rational systems): Let σ : R → R be a globally Lipschitz function which satisfies (A1) and assume that U is finite. Consider an input-output map p : U pc → P C([0; +∞[, R p ). If Σ is a RNN with activation function σ and input space U, and Σ is a realization of p, then there exists a rational system which is a realization of p. The proof of Theorem 1 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF]. The proof of the theorem relies on defining a rational system associated with the RNN. In order to avoid excessive notation, in the sequel we identify the sum of fraction of multi-variable polynomials

N k=1

Pi Qi with the fraction obtained by bringing all summands to the same denominator, i.e., N k=1 PiΠ N r=1,r =k Qi Π N k=1 Qi . Definition 7: Let Σ be an RNN, whose activation function satisfies (A1), and assume that ξ = (ξ 1 , . . . , ξ N ) T satisfies (4). Define the rational system R(Σ) associated with the RNN Σ = (A, B, C, U, σ, x 0 ), A = (a i,j ) n i,j=1 , C = (c k,i ) k=1,...,p,i=1,...,n as follows: The proof of Theorem 1 relies on the following simple result, which is interesting on its own right. Let ξ = (ξ 1 , . . . , ξ n ) is from (4), and define F : R n → R L , as F (x) = (z 1 , . . . , z nN |U | , x T ) T , where z φ(i,j,α) = ξ i (e T j (Ax + Bα)) and for any i = 1, . . . , N, j = 1, . . . , n, α ∈ U, φ(i, j, α)

∀i = 1, . . . , N, j = 1, . . . , n, α ∈ U : υi,j,α (t) = U i (υ j,α (t)) V i (υ j,α (t)) { n l=1 a j,l U 0 (υ l,β (t)) V 0 (υ l,β (t)) } if u(t) = β υ j,α (t) = (υ 1,j,α (t), . . . , υ N,j,α (t)), υ j,α (0) = ξ e T j (Ax 0 + Bα)) ẋj (t) = U 0 (υ j,β (t)) V 0 (υ j,β ) if u(t) = β,
= N • |U | • (j -1) + N • (r -1) + i if α = α r , r = 1, . . . , K.
Lemma 2: If (x, u, y) is a solution of the RNN Σ = (A, B, C, U, σ, x 0 ), then (υ, u, y), with υ(t) = F (x(t)) for t ≥ 0, is a solution of R(Σ). The proof of Lemma 2 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

IV. EXISTENCE OF AN RNN REALIZATION

Theorem 1 allows us to formulate a necessary condition for realizability of an input-output map by an RNN, using conditions of [START_REF] Němcová | Realization theory for rational systems: The existence of rational realizations[END_REF]Theorem 5.16] for existence of a realization by a rational system. In order to present this condition, we need to introduce additional notation and terminology.

Definition 8: An input-output map p :

U pc → P C([0; +∞[, R p ) is causal if, for all t
0 and for all u, v ∈ U pc such that ∀s ∈ [0, t] : u(s) = v(s), it holds that ∀s ∈ [0, t] : p(u)(s) = p(v)(s).

In other words, p is causal, if p(u)(t) depends only on the values of u on the interval [0, t]. If p is the input-output map of a control system, then causality must necessarily hold. Definition 9: An input-output map p : U pc → P C([0; +∞[, R p ) is analytic if, for all k ∈ {1, . . . , p} and for all α 1 , . . . , α l ∈ U, l > 0, the function φ p,k,α1,...,α l : ([0, +∞[) l → R is analytic, where φ p,k,α1,...,α l (t 1 , . . . , t l ) = p k (u α1,...,α l t1,...,t l )(T l ) u α1,...,α l t1,...,t l (t) =

α i if t ∈ [T i-1 , T i [, i = 1, . . . , l α l if t ≥ T l T 0 = 0 , T i = i j=1 t j , i = 1, 2, . . . , l. (6)
where p k is the k-th component of p, i.e. p(u) = (p 1 (u), . . . , p p (u)) T for all u ∈ U pc . Notation 3: We denote by A(U pc ) the set of causal analytic input-output maps of the form p : U pc → P C([0; +∞[, R). Note that the set A(U pc ) forms an algebra over the field of real numbers with the usual point-wise addition, multiplication and multiplication by scalar. It can be shown that the algebra A(U pc ) is isomorphic to the ring of functions A U pc → R defined in [START_REF] Němcová | Realization theory for rational systems: The existence of rational realizations[END_REF]Definition 4.3], therefore it is an integral domain.

Definition 10: Let ϕ ∈ A(U pc ) and define the derivative D α ϕ of ϕ along α ∈ U as the function D α ϕ : U pc → P C([0; +∞[, R), such that for all u ∈ U pc , for all t 0,

D α ϕ(u) (t) = d ds ϕ u α )(t + s) | s=0 , u α (τ ) = u(τ ) τ ∈ [0, t[ α τ > t
It is easy to see that D α ϕ is also causal and analytic, and hence D α ϕ belongs to A(U pc ). Now we define the observation algebra of an input-output map.

Definition 11: Let p : U pc → P C([0; +∞[, R p ) be an analytic and causal input-output map. The observation algebra of p, denoted by A obs (p), is the smallest sub-algebra of the algebra A U pc such that the following holds.

• Consider the components p k :

U pc → P C([0; +∞[, R), k = 1, . . . , p of p, i.e., ∀u ∈ U pc : p(u) = (p 1 (u), . . . , p p (u)) T . For every k = 1, . . . , p, p k ∈ A obs (p). • For every g ∈ A obs (p), D α g ∈ A obs (p), α ∈ U, i.e.,
A obs (p) is closed under taking derivatives D α , α ∈ U. We call the observation field, denoted by Q obs (p), the field of fractions of A obs (p).

If an input-output map p : U pc → P C([0; +∞[, R p ) is realized by an RNN Σ = (A, B, C, U, σ, x 0 ), with σ satisfying Assumption (A1), then the rational system R(Σ), given in Definition 7, also realizes the input-output map p by Theorem 1 and Lemma 1. Thus the transcendence degree of the observation algebra A obs (p) is finite by [START_REF] Němcová | Realization theory for rational systems: The existence of rational realizations[END_REF]Theorem 5.16]. We are now in the position to state a necessary condition for existence of an RNN realization, which summarizes the arguments above.

Theorem 2 (Existence of an RNN: necessary condition): The input-output map p : U pc → P C([0; +∞[, R p ) has a realization by an RNN whose activation function satisfies Assumption (A1) only if p is causal, analytic and trdeg A obs (p) < +∞.

V. MINIMALITY, REACHABILITY AND OBSERVABILITY

In this section, we first provide sufficient conditions for minimality of RNNs. We also provide a Hankel-rank like condition for minimality of the RNNs. Finally we relate reachability and observability properties for RNNs with the corresponding properties of rational systems.

A. Sufficient conditions for minimality of RNNs

Define the dimension of a RNN Σ of the form (1) (respectively rational system R of the form (2)) as the number n of its state variables, and we denote the dimension of Σ (respectively R) by dim(Σ) (respectively dim(R)) . Note that our definition of dimension of a rational system coincides with that of [START_REF] Němcová | Realization theory for rational systems: Minimal rational realizations[END_REF]Definition 13]. We say that RNN Σ with activation function σ is a σ-minimal realization (respectively the rational system R is a minimal realization) of an input-output map p, if Σ (respectively R) is a realization of p and there exists no RNN Σ with activation function σ (respectively no rational system R ), such that Σ (respectively R ) is a realization of p and dim(Σ ) < dim(Σ) (respectively dim(R ) < dim(R)). . Lemma 3: Let Σ = (A, B, C, σ, x 0 ) be an RNN, whose activation function σ satisfies (A1) and assume that Σ is a realization of the input-output map p. If the rational system R(Σ), given in Definition 7, is a minimal realization of p, then Σ is a minimal RNN realization of p. The proof of Lemma 3 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

Unfortunately, in most of the cases, R(Σ) will not be minimal. Intuitively, this has to do with the fact that the states x 1 (t), . . . , x n (t) of R(Σ) are integrals of the other states, leading to lack of observability if n > 1 and x(0) is chosen so that Cx(0) = 0. In order to remedy this problem, we define auxiliary rational systems associated with RNNs.

Definition 12: Let Σ = (A, B, C, U, σ, x 0 ) be an RNN, whose activation function satisfies (A1), assume that ξ = (ξ 1 , . . . , ξ N ) satisfies (4). Define the auxiliary rational system R aux (Σ) associated with the RNN Σ, A = (a i,j ) n i,j=1 , C = (c k,i ) k=1,...,p,i=1,...,n as follows: ∀i = 1, . . . , N, j = 1, . . . , n, α ∈ U : υi,j,α (t) = U i (υ j,α (t))

V i (υ j,α (t)) { n l=1 a j,l U 0 (υ l,β (t)) V 0 (υ l,β (t)) } , if u(t) = β υ j,α (t) = (υ 1,j,α (t), . . . , υ n,j,α (t)) T , υ j,α (0) = ξ e T j (Ax 0 + Bα) , y k,α (t) = n i=1 c k,i U 0 (υ i,α (t)) V 0 (υ i,α (t)) , k = 1, . . . , p Remark 3 (Polynomial R aux (Σ)): If Assumption (A1
) is satisfied with polynomial equations, i.e. V i = 1, i = 0, . . . , N , like in examples from Example 1, then R aux (Σ) is a polynomial system.

Note that the rational system R aux (Σ) does not realize the input-output map p Σ,x0 of the RNN Σ, but it realizes an input-output map which is constructed from derivatives of p .More precisely, let p : U If the RNN Σ realizes p, then the rational system R aux (Σ) realizes the input-output map p defined in [START_REF] Dalessandro | Realization and structure theory of bilinear dynamical systems[END_REF]. The proof of Lemma 4 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

Lemma 5: Let Σ = (A, B, C, U, σ, x 0 ) be an RNN and let p be an input-output map realized by Σ. If R aux (Σ) is a minimal realization of p, then the RNN Σ is σ-minimal realization of p. The proof of Lemma 5 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

It is known from [START_REF] Němcová | Rational Systems in Control and System Theory[END_REF], [START_REF] Němcová | Realization theory of Nash systems[END_REF] that algebraic, rational and semi-algebraic observability and algebraic reachability characterize minimality of rational systems and these properties can be checked using methods of computational algebra [START_REF] Němcová | Observability reduction algorithm for rational systems[END_REF]. In particular, we can derive sufficient conditions for the minimality of the RNN Σ using these reachability and observability concepts for rational systems. To this end, below the notions of algebraic reachability and algebraic/semialgebraic observability for rational systems. Let R be a rational system of the form (2). Define the set R R (υ 0 ) of reachable states of R as the set formed by υ(t) such that t 0 and (υ, u, y) is a solution of R with υ(0) = υ 0 . The system R is said to be algebraically reachable, if there is no non-trivial polynomial which is zero on R R (υ 0 ). The system R is called accessible, if R R (υ 0 ) contains an open subset of R n . It is clear that accessibility of implies algebraic reachability.

Let R be a rational system of the form (2). Recall from [START_REF] Němcová | Rational Systems in Control and System Theory[END_REF]Definition 3.19] or from [START_REF] Bartoszewicz | Rational systems and observation fields[END_REF]Definition 4] that observation algebra of R, denoted by A obs (R), is the smallest sub-algebra of the field of rational functions R(X 1 , . . . , X n ) which contains h k,1 h k,2 , k = 1, . . . , p and which is closed under taking the formal Lie derivatives with respect to the formal vector fields

f α = n i=1 Pi,α Qi,α ∂ ∂Xi . If R is polynomial, i.e. Q i,α = 1, i = 1, . . . , n, h k,2 = 1, k = 1, . . . , p, then A obs (R) is the sub-algebra of the ring of polynomials R[X 1 , . . . , X n ].
Following [START_REF] Němcová | Rational Systems in Control and System Theory[END_REF] we say that the rational system R is algebraically observable, if A obs (R) = R[X 1 , . . . , X n ]. Following [START_REF] Němcová | Realization theory of Nash systems[END_REF] that R is semi-algebraically observable if trdeg(A obs (R)) = n. We will say that R is observable, if for every two distinct initial states υ " 0 , υ 0 there exists solutions (υ " , u, y " ) and (υ , u, y ) of R such that υ " (0) = υ " 0 , υ (0) = υ 0 , and y " = y . It is easy to see that algebraic observability implies semi-algebraic observability. Moreover, for polynomial systems algebraic observability implies observability [START_REF] Bartoszewicz | Minimal polynomial realizations[END_REF].

Lemma 6 (Sufficient conditions): If one of the conditions below holds, then Σ is σ-minimal realization of p:

• R aux (Σ) is semi-algebraically observable and algebraically reachable. • R aux (Σ) is polynomial, it is algebraically observable and algebraically reachable. • R aux (Σ) is polynomial, it is algebraically observable and accessible. The above lemma is then a direct consequence of Lemma 5 and [20,Proposition 6], its detailed proof is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

B. A Hankel-rank like condition for minimality of RNNs

Recall that a linear system is a minimal realization of its input-output map, if and only if the dimension of this system equals the rank of the Hankel-matrix of the inputoutput map. We would like to formulate a similar result, where the role of the rank of the Hankel-matrix is played by the observation algebra A obs (p). To this end, recall from [20, Lemma 1, Theorem 4] that R(Σ) is minimal if and only if dim(R(Σ)) = trdegA obs (p). In a similar manner, if A obs (p) is the observation algebra of the input-output map p defined in [START_REF] Dalessandro | Realization and structure theory of bilinear dynamical systems[END_REF], then R aux (Σ) is minimal if and only if dim(R aux (Σ)) = trdegA obs (p).

Lemma 7: Assume that σ satisfies (A1), and let the RNN Σ be a realization of the input-output map p. If one of the following conditions hold for n = dim(Σ): (1) trdegA obs (p) = n(1 + |U|N ), or (2) trdegA obs (p) = n|U|N , then Σ is a σminimal realization of p. The proof of Lemma 7 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

C. Some aspects of reachability and observability of RNNs

One may wonder how restrictive the conditions of Lemma 6 are, and how they relate to accessibility/reachability and observability of the RNN Σ studied in [START_REF] Albertini | State observability in recurrent neural networks[END_REF], [START_REF] Qiao | Further results on controllability of recurrent neural networks[END_REF], [START_REF] Albertini | Forward accessibility for recurrent neural networks[END_REF]. In order to present this relationship more precisely, we introduce the following terminology. Define the reachable set of an RNN Σ = (A, B, C, U, σ, x 0 ) R Σ (x 0 ) = {x(t) | t 0 , (x, u, y) is a solution of Σ, x(0) = x 0 }. We will say that Σ is accessible, if R Σ (x 0 ) contains an open subset of R n , we say that Σ is algebraically reachable if there is no non-trivial polynomial which is zero on R Σ (x 0 ). We say that Σ is span-reachable, if the linear span of the elements R Σ (x 0 ) is R n , i.e. Σ is reachable if there exist no linear function which is zero on R Σ (x 0 ). Clearly, if Σ is accessible, then it is algebraically reachable, and if Σ is algebraically reachable, then it is span-reachable. We say that the RNN Σ is weakly observable if for every initial state x ∈ R n there is an open subset V of R n such that x ∈ V and for every x = x ∈ V , there exist solution (x, u, y) and (x , u, y ) of Σ, with x(0) = x and x (0) = x, such that

Remark 2 :

 2 and x j (0) = e T j x 0 ,y k (t) = n i=1 c k,i x i (t), k = 1, . . . , p If Assumption (A1) is satisfied with polynomial equations, i.e. V i = 1, i = 0, . . . , N , like in examples from Example 1, then R(Σ) is a polynomial system.

  pc → P C([0; +∞[, R p ) be an inputoutput map realized by the RNN Σ. Define the input-output map p : U pc → P C([0; +∞[, R p|U | ) as follows:∀u ∈ U pc , p(u) = (D α1 p(u), . . . , D α K p(u)) T ,(7)where,D α p(u) = (D α p 1 (u), . . . , D α p p (u)) T , α ∈ U, D α p k (u)is defined in Definition 10 and p(u) = (p 1 (u), . . . , p p (u)) T . Lemma 4: Let Σ = (A, B, C, U, σ, x 0 ) be an RNN and let p : U pc → P C([0; +∞[, R p ) be an input-output map.

VI. CONCLUSIONS

We have shown that input-output maps of a large class of recurrent neural networks can be represented by rational/polynomial systems, and we used this fact to derive necessary and sufficient conditions for existence of a realization by a recurrent neural network and its minimality. Future research will be directed towards deriving a more complete realization theory of recurrent neural network and for using the results of realization theory for analyzing machine learning algorithms.

y = y . Then the RNN Σ is observable if for every initial state x ∈ R n , V = R n in the latter definition.

Lemma 8: Let Σ = (A, B, C, U, σ, x 0 ) be an RNN.

• If R aux (Σ) is algebraically reachable, then Σ is spanreachable. In particular, if R aux (Σ) is accessible, then Σ is span-reachable. [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF]. Note that invertability of σ holds for many commonly used activation functions, see Example 1. Observe that accessibility, and algebraic / semi-algebraically observability conditions for rational / polynomial systems can be checked by using methods of computer algebra [START_REF] Němcová | Observability reduction algorithm for rational systems[END_REF].

From [START_REF] Jakubczyk | Realization theory for nonlinear systems: three approaches[END_REF]Theorem 1.12] it follows that if the RNN Σ is accessible and weakly observable, then it is a minimal realization of its input-output map p. From the comparison between the conditions of Lemma 6 with those of Lemma 8 it is clear that minimality of R aux (Σ) is a much weaker condition than accessibility and weak observability of Σ. This suggests that using realization theory of rational systems is likely to yield more useful results for RNNs than using realization theory of general analytic systems.

Recall from [3, Theorem 1] that a necessary condition for observability of Σ = (A, B, C, U, σ, x 0 ) is that the largest A-invariant coordinate subspace of Σ included in Ker(C) is trivial. More precisely, following [START_REF] Albertini | State observability in recurrent neural networks[END_REF] we say that a vector subspace V of R n is a coordinate subspace if it is spanned by some vectors from the canonical basis of R n , i.e. there exists an integer s > 0 and integers i 1 , . . . , i s ∈ {1, . . . , n} such that V is spanned by e i1 , . . . , e is , where (e 1 , . . . , e n ) denotes the canonical basis of R n . Let O c (A, C) be the largest coordinate subspace which is A-invariant and contained in Ker(C).

Lemma 9: If R aux (Σ) is polynomial and it is semialgebraically observable, then there exists no non-trivial coordinate subspace which is A-invariant and contained in Ker(C), i.e. O c (A, C) = {0}. The proof of Lemma 9 is presented in [START_REF] Defourneau | Realization theory of recurrent neural networks and rational systems[END_REF].

By [3, Theorem 1] O c (A, C) = {0}, it is also sufficient if ker(C) ∩ ker(A) = {0}, the activation function σ satisfies only the IPP property, given in [START_REF] Albertini | State observability in recurrent neural networks[END_REF] for example, and if B verifies a condition on its rows. But here we do not need the latter hypothesis on Σ.


Manas Mejari
email: manas.mejari@idsia.ch

Mihály Petreczky
email: mihaly.petreczky@ec-lille.fr

Realization and identification algorithm for stochastic LPV state-space models with exogenous inputs

In this paper, we present a realization and an identification algorithm for stochastic Linear Parameter-Varying State-Space Affine (LPV-SSA) representations. The proposed realization algorithm combines the deterministic LPV input output to LPV state-space realization scheme based on correlation analysis with a stochastic covariance realization algorithm. Based on this realization algorithm, a computationally efficient and statistically consistent identification algorithm is proposed to estimate the LPV model matrices, which are computed from the empirical covariance matrices of outputs, inputs and scheduling signal observations. The effectiveness of the proposed algorithm is shown via a numerical case study. 71 pt 0.986 in 25 mm 43 pt 0.597 in 15.2 mm 43 pt 0.597 in 15.2 mm 71 pt 0.986 in 25 mm 71 pt 0.986 in 25 mm 43 pt 0.597 in 15.2 mm 43 pt 0.597 in 15.2 mm 71 pt 0.986 in 25 mm

INTRODUCTION

Identification of Linear Parameter-Varying (LPV) models has gained significant attention over the past few years, owing to their ability to describe the behavior of many time-varying and non-linear systems. Many approaches have been proposed for the identification of LPV models, in input-output [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] Mejari | A bias-correction method for closed-loop identification of Linear Parameter-Varying systems[END_REF][START_REF] Piga | LPV system identification under noise corrupted scheduling and output signal observations[END_REF] as well as State-Space (SS) representations [START_REF] Felici | Subspace identification of MIMO LPV systems using a periodic scheduling sequence[END_REF][START_REF] Tanelli | Identification of LPV state space models for autonomic web service systems[END_REF][START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF][START_REF] Verdult | Kernel methods for subspace identification of multivariable LPV and bilinear systems[END_REF]. The reader is referred to [START_REF] Tóth | Modeling and identification of linear parameter-varying systems[END_REF] for a detailed summary of the available LPV identification approaches.

Controller design approaches often require the LPV models to be in SS representation with an affine dependency on the scheduling variable. To this end, realization theory of LPV models plays a key role in understanding the conditions under which the observed behavior of a system can be realized by a state-space affine representation. It also allows to formulate identification algorithms for estimating state-space representation from a finite set of observations. The realization theory for deterministic Linear Parameter-Varying State-Space with Affine dependence (LPV-SSA) representation has been developed in [START_REF] Tóth | On the state-space realization of LPV input-output models: Practical approaches[END_REF]; [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF]. The results of [START_REF] Tóth | On the state-space realization of LPV input-output models: Practical approaches[END_REF]; [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF] were used to derive LPV-SS identification algorithm in [START_REF] Cox | Estimation of LPV-SS models with static dependency using correlation analysis[END_REF][START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF]. These methods are focused on deterministic realizations, which for certain control and filtering problems are too restrictive. In this paper, This work was partially funded by CPER Data project, co-financed by European Union with the financial support of European Regional Development Fund, French State and the French Region of Hauts-de-France, and by CNRS project PEPS blanc BayesRealForRNN.

we focus on formulating a realization algorithm and a related identification algorithm for stochastic LPV-SSA representations. The main idea is to decompose the stochastic LPV-SSA realization/identification problem into two independent problems: realization/identification of deterministic part which depends only on the input, and realization/identification of stochastic part. To this end, the proposed algorithm is based on the combination of correlation analysis [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF] for deterministic realization and stochastic covariance identification algorithm for stochastic LPV-SSA representations (Mejari and Petreczky, 2019a).

The algorithm presented in this paper extends the results of Mejari and Petreczky (2019a), to the case of stochastic LPV-SSA representations with exogenous inputs. The proposed approach differs significantly from the subspace based identification methods for stochastic LPV-SSA representations (van Wingerden and Verhaegen, 2009;[START_REF] Dos Santos | Identification of bilinear systems with white noise inputs: An iterative deterministic-stochastic subspace approach[END_REF][START_REF] Favoreel | Subspace identification of bilinear systems subject to white inputs[END_REF]. First, the cited papers do not deal with the realization problem. In particular, while the possibility of decomposing the output into a deterministic and purely stochastic components is sometimes claimed in the literature, the formal details of such a decomposition were never addressed. Second, in contrast to the literature mentioned above, the identification algorithm proposed in this paper is provenly consistent and it does not require local observability assumptions. The downside is that the proposed algorithm is provenly consistent only for a specific class of scheduling signals and stochastic LPV-SSA representations. Moreover, the proposed algorithm avoids the curse of dimensionality, but this comes at a price of either using some prior knowledge on the system to determine the correct selection of the rows and columns of a Hankel-matrix or using an exhaustive search to find such a selection.

Margin requirements for the other pages Paper size this page A4

The paper is organized as follows. In Section 2, we present the problem formulation. Section 3 presents the formal definition and basic properties of the class of LPV state-space representations considered in this paper. In Section 4, we formalize the decomposition of outputs of such LPV state-space representations into stochastic and deterministic components. In Section 5, we present the realization algorithm for stochastic LPV state-space representations, and in Section 6 we present the related identification algorithm. Finally, in Section 7 we illustrate the results with a numerical example.

PROBLEM FORMULATION

Let y, u, µ be stochastic processes taking values in R ny , R nu and R nµ , representing the output process, the input process, and the scheduling signal process respectively. We define a discrete-time Linear Parameter-Varying State-Space Affine (LPV-SSA) representation of the process (y, u, µ) as the discrete-time system of the form

x(t + 1) = nµ i=1 (A i x(t) + B i u(t) + K i v(t))µ i (t), y(t) = Cx(t) + Du(t) + v(t),
(1) where,

A i ∈ R nx×nx , B i ∈ R nx×nu , K i ∈ R nx×ny , ∀i = 1, . . . , n µ , C ∈ R ny×nx and D ∈ R ny×nu are real constant matrices, and v is a white noise process, i.e., E[v(t)v T (s)] = 0, s = t and E[v(t)v T (t)µ i (t)] = Q i > 0, i = 1, . . . , n µ .
The realization and identification problems considered in this paper are as follows. Problem 1. (Realization problem). For process (y, u, µ), find matrices ({A i , B i , K i } nµ i=1 , C, D) and processes x, v such that (1) is a representation of (y, u, µ). Problem 2. (Identification problem). Assume that y : Z → R ny is a sample path of the output process y, u : Z → R nu is a sample path of the input process u and µ : Z → R nµ is a sample path of the scheduling process µ, corresponding to the same random event ω ∈ Ω. Given a dataset {y(t), u(t), µ(t)} N t=1 consisting of N samples of the output, input and scheduling process, compute from this dataset the estimates

{{ ÂN i , BN i , KN i , QN i } nµ i=1 , ĈN , DN }, such that as N → ∞, the estimated matrices {{ ÂN i , BN i , KN i , QN i } nµ i=1 , ĈN , DN } converge to matrices {{A i , B i , K i , Q i } nµ i=1 , C, D} such that the LPV-SSA (1) with Q i = E[v(t)v (t)µ 2 i (t)], i = 1, . . . , n µ
, is a representation of (y, u, µ).

PROPERTIES OF LPV-SSA REPRESENTATION

In order to make Problems 1-2 well-posed, we have to impose additional constraints on the class of processes (y, u, µ) and on the class of LPV-SSA representations. To this end, we need the following notation and terminology. Notation 1. (Σ). Let Σ = {1, . . . , n µ }.

A non empty word over Σ is a finite sequence of letters, i.e.,

w = σ 1 σ 2 • • • σ k , where 0 < k ∈ Z, σ 1 , σ 2 , . . . , σ k ∈ Σ.
The set of all nonempty words is denoted by Σ + . We denote an empty word by . Let Σ * = ∪ Σ + . The concatenation of two nonempty words

v = a 1 a 2 • • • a m and w = b 1 b 2 • • • b n is defined as vw = a 1 • • • a m b 1 • • • b n for some m, n > 0. Note that if w = or v = , then v = v and w = w, moreover, = .
The length of the word w ∈ Σ * is denoted by |w|, and | | = 0. Assumption 1. (White noise scheduling). The scheduling process µ = [1, µ 2 , . . . , µ nµ ] T is zero-mean independent identically distributed (i.i.d.) such that, for all t ∈ Z, we have µ 1 (t) ≡ 1, and for each σ = 2, . . . , n µ , µ σ is a zero mean i.i.d. process.

We define scalars E[µ 2 σ (t)] = p σ , for all t ∈ Z. In particular, p 1 = 1. For every word w ∈ Σ + where w

= σ 1 σ 2 • • • σ k , k ≥ 1, σ 1 , . . . , σ k ∈ Σ,
we define the process µ w and the number p w as follows

µ w (t) = µ σ1 (t -k + 1)µ σ2 (t -k + 2) • • • µ σ k (t), ∀t ∈ Z p w = p σ1 p σ2 • • • p σ k .
We set µ (t) = 1 and p = 1. For a process r ∈ R nu , for each w ∈ Σ + we define the process z r w as

z r w (t) = r(t -|w|)µ w (t -1) 1 √ p w , ∀t ∈ Z, (2)
which is interpreted as the past of r w.r.t. {µ σ } σ∈Σ .

We remark that in this paper, the notion of Zero Mean Wide Sense Stationary w.r.t. Inputs (ZMWSSI) process and Square Integrable w.r.t. Input (SII) process, will be a central notion for the mathematical framework of stochastic LPV-SSA representations. We refer the reader to [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF] for the detailed definition of ZMWSSI and SII processes.

All the process considered in this paper will be assumed to be ZMWSSI and SII process w.r.

(t)(z r v (t)) T ] = 0, v = w, E[z r σw (t)(z r σw (t)) T ] = E[z r σ (t)(z r σ (t)) T ] > 0, for all w, v ∈ Σ + .
Using the concept of ZMWSSI process and white noise process w.r.t. µ, we can formulate the main assumption regarding the processes (y, u, µ). Assumption 2. Assume that µ satisfies Assumption 1, and y T u T T is a ZMWSSI and SII process w.r.t. µ, and u is a white noise process w.r.t. µ, and the covariance

E[z u σ (t)(z u σ (t)) T ] = E[u(t -1)(u(t -1)) T ] = Λ u > 0 does not depend on σ ∈ Σ.
Next, we recall from Mejari and Petreczky (2019a) the notion of a stationary stochastic LPV-SSA representation of a process r without inputs. Definition 2. A stationary LPV-SSA representation without inputs of a process r taking values in

R p , is a tuple ({ Ãσ , Kσ } nµ σ=1 , C, D, x, ṽ), where Ãσ ∈ R ñ×ñ , Kσ ∈ R ñ× m, C ∈ R p×ñ and v is a process taking values in R m such that such that 1. xT ṽT T is a ZMWSSI process, and E[x(t)(z ṽ w (t)
) T ] = 0 for all w ∈ Σ + . 2. ṽ is a white noise process w.r.t. µ.

The eigenvalues of the matrix

σ∈Σ p σ Ãσ ⊗ Ãσ are inside the open unit circle. 4. x(t+1) = nµ i=1 ( Ãi x(t)+ Ki ṽ(t))µ i (t), r(t) = C x(t)+ Dṽ(t).
We call x the state process and ṽ the noise process. In the terminology of [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF], a stationary LPV-SSA without inputs u, corresponds to a stationary generalized bilinear system w.r.t. the scheduling inputs {µ σ } σ∈Σ . If a process r has a stationary LPV-SSA representation without inputs, then r is a ZMWSSI process and x is uniquely determined by ṽ and the matrices ( C, D, { Ãσ , Kσ } σ∈Σ ). In order to define this notion more precisely, let us introduce the following notation. Notation 2. (Matrix Product). Consider a collection of square matrices

A σ ∈ R n×n , σ ∈ Σ. For any word w ∈ Σ + of the form w = σ 1 σ 2 • • • σ k , k > 0 and σ 1 , . . . , σ k ∈ Σ, we define A w =A σ k • • • A σ2 A σ1 . For an empty word , A = I n .
From [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF] it follows that

x(t) = σ∈Σ,w∈Σ * √ p σw Ãw Kσ z ṽ σw (t), (3)
where the infinite sum on the right-hand side is absolutely convergent in the mean square sense.

Using the notion of a stationary LPV-SSA without inputs, we can define the class of LPV-SSA representation with inputs which will be considered in this paper.

Definition 3. (Stationary LPV-SSA). The LPV-SSA representation (1) is stationary with input u, if ({A σ , [K σ B σ ]} σ∈Σ , C, I ny D x, v T u T T )
is a stationary LPV-SSA representation of y without inputs as in Definition 2, and the orthogonality condition

E[v(t)u T (t)µ 2 σ (t)] = 0, ∀σ ∈ Σ holds.
From (3) it follows that for a stationary LPV-SSA representation with input u of the form (1),

x(t) = w∈Σ * ,σ∈Σ √ p σw A w (K σ z v σw (t) + B σ z u σw (t)) ,
where the infinite sums on the right hand side are absolutely convergent in the mean-square sense. That is, the matrices and the noise processes determine the state process of a stationary LPV-SSA (with or without inputs) uniquely.

DECOMPOSITION OF THE OUTPUT OF LPV-SSA REPRESENTATION

It turns out that the output process of stationary LPV-SSA representations admits a decomposition into deterministic and stochastic parts. The deterministic part depends only on the input process, while the stochastic part depends only on the noise process. This decomposition does not depend on the particular choice of LPV-SSA representation, but only on the output process at hand. Definition 4. (Deterministic and stochastic components).

Assume the processes (y, u, µ) satisfy Assumption 2. Define the deterministic component y d of y as follows

y d (t) = E l [y | {z u w (t)} w∈Σ + ∪ {u(t)}], (4
) where E l [•] is the orthogonal projection as defined in [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]. Define the stochastic component of y as y s (t) = y(t) -y d (t).

(5)

From the definition it follows that y(t) = y d (t) + y s (t),

i.e., the process y(t) can be represented as the sum of its deterministic and stochastic components. In case when the process admits an LPV-SSA representation, the stochastic and deterministic components satisfy the following properties. Lemma 1. (Decomposition of y). Assume that there exists a stationary LPV-SSA representation of (y, u, µ) of the form (1) and that (y, u, µ) satisfy Assumption 2. It then follows that

x d (t + 1) = nµ i=1 (A i x d (t)+B i u(t))µ i (t), y d (t) = Cx d (t)+Du(t), (6)
and ({A σ , B σ } σ∈Σ , C, D, x d , u) is a stationary LPV-SSA representation of y d without inputs and with noise process u, moreover,

x s (t + 1) = nµ i=1 (A i x s (t)+K i v(t))µ i (t),
y s (t) = Cx s (t)+v(t), ( 7) and ({A σ , K σ } σ∈Σ , C, I ny , x s , v) is a stationary LPV-SSA representation of y s without inputs, where

x d (t) = E l [x(t) | {z u w (t)} w∈Σ + ∪ {u(t)}] (8) x s (t) = x(t) -x d (t) (9)
The proof of Lemma 1 is presented in [START_REF] Mejari | Realization and identification algorithm for stochastic LPV state-space models with exogenous inputs[END_REF]. Thus, y s depends only on the noise v, and y d does not depend on the noise but it depends only on input u.

In fact, the converse of Lemma 1 also holds. ) is a stationary LPV-SSA representation of y s without inputs in forward innovation form, i.e., assume that the process e s is the so called innovation process of y s as defined in [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]:

e s (t) = y s (t) -E l [y s (t) | {z y s w (t)} w∈Σ + ] (10)
Then, tuple ({ Âi , Ki , Bi } nµ i=1 , Ĉ, D, x, e s ) is a stationary LPV-SSA representation of y with input u, where

x(t) = (x d (t)) T (x s (t)) T T Âσ = diag( Âd σ , Âs σ ), Bσ = ( Bd σ ) T 0 T nx×nu T Kσ = 0 T nx×ny ( Ks σ ) T T , Ĉ = Ĉd Ĉs , D = Dd . (11)
Moreover, the innovation process e s satisfies e s (t) = y(t)

-E l [y(t) | {z y w (t), z u w (t)} w∈Σ + ∪ {u(t)}] (12)
The proof of Lemma 2 is presented in [START_REF] Mejari | Realization and identification algorithm for stochastic LPV state-space models with exogenous inputs[END_REF]. Thus, the problem of realization of y can be decomposed into two problems:

P1 finding a stationary LPV-SSA representation Σ d without inputs of y d , such that the noise process of Σ d is u, P2 finding a stationary LPV-SSA representation Σ s without inputs of y s = yy d , such that the noise Margin requirements for the other pages Paper size this page A4 process e s of Σ s is the innovation process of y s as defined in Mejari and Petreczky (2019a).

Moreover, the innovation process e s (t) is the error of projecting y(t) onto the linear space spanned by the products of the past values of y, u and the scheduling process µ, as defined in (12).

REALIZATION ALGORITHMS

In this section, we first recall the basis reduced Ho-Kalman realization algorithm for deterministic LPV state-space representations. In turn, this algorithm will be used for covariance realization algorithms for estimating LPV-SSA representations of y d , y s , presented in Section 5.2-5.3.

Basis reduced Ho-Kalman realization algorithm

Recall from [START_REF] Petreczky | Realization theory for LPV state-space representations with affine dependence[END_REF]; [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF] that a deterministic LPV-SSA representation (with affine dependence) is a system of the form

x(t + 1) = nµ i=1 (A i x(t) + B i u(t))µ i (t), y(t) = Cx(t) + Du(t), (13)
where A i , B i , C, D are matrices of suitable dimensions, x : Z → R nx is the state trajectory u : Z → R nu is the input trajectory y : Z → R ny is the output trajectory. In order to avoid technical problems, we assume that x, u, y all have finite support, i.e. there exist a t 0 ∈ Z, such that x(s) = 0, y(s) = 0, u(s) = 0 for all s < t 0 . We identify a deterministic LPV-SSA of the form ( 13

. α ⊆ Σ n × {1, 2, • • • , n y } and β ⊆ Σ × Σ n × {1, 2, • • • , n u } 2. card(α) = card(β) = n,
where card denotes cardinality of the set.

We will fix the following ordering of α and β.

α = {(u i , k i )} n i=1 , β = {(σ j , v j , l j )} n j=1 , (15)
u i ∈ Σ n , k i ∈ {1, 2, • • • , n y }, σ j ∈ Σ, v j ∈ Σ n , l j ∈ {1, 2, • • • , n u } Example 1.
Consider n = 2, number of outputs and inputs n y = n u = = 2, and scheduling signal dimension n µ = 2, we have, Σ n = { , 1, 2, 11, 12, 21, 22}. Then, one of the n-selection pair (α, β) can be chosen as, for e.g.,

α = {(u 1 , k 1 ) , (u 2 , k 2 )} = {( , 1) , (11, 2)} and β = {(σ 1 , v 1 , l 1 ) , (σ 2 , v 2 , l 2 )} = {(1, 21, 1) , (2, 22, 2)}.
Let M : Σ * → R ny×nu be a map, values of which represent potential sub-Markov parameters ( 14) of an LPV-SSA. Let us now define the Hankel matrix H M α,β ∈ R n×n as follows: i, j = 1, . . . , n, the (i, j)-th element of H M α,β is of the form

H M α,β i,j = [M (σ j v j u i )] ki,lj , (16)
[M (σ j v j u i )] ki,lj denotes the entry of M (σ j v j u i ) on the k i -th row and l j -th column, and (u i , k i ) ∈ α, (σ j , v j , l j ) ∈ β are as in the ordering of ( 15).

In addition, we define the σ-shifted Hankel-matrix H M σ,α,β ∈ R n×n as follows: its i, j-th entry is given by

H M σ,α,β i,j = [M (σ j v j σu i )] ki,lj . (17)
Moreover, let us define Hankel matrices

H M α,σ ∈ R n×nu and H M β ∈ R ny×n as follows H M α,σ i,j = [M (σu i )] ki,j , j = 1, . . . , n u (18)
H M β i,j = [M (σ j v j )] i,lj , i = 1, . . . , n y (19)
Consider the model matrix computations summarized in Algorithm 1, using Hankel matrices and selections.

= (H M α,β ) -1 H M σ,α,β , Bσ = (H M α,β ) -1 H M α,σ , Ĉ = H M β 2. Compute D = M ( ).
Output: Matrices ({ Âσ , Bσ } σ∈Σ , Ĉ, D)

Lemma 3. (Adapted from [START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF]). Let the (n, n y , n u )-selection (α, β) be such that rank(H M α,β ) = n, and assume that there exists a deterministic LPV-SSA representation S * of dimension n such that M = M S * . Then the tuple Ŝ = ({ Âσ , Bσ } σ∈Σ , Ĉ, D), returned by Algorithm 1, when applied to the matrices 16)-( 19)) and M ( ), is a minimal dimensional deterministic LPV-SSA representation such that M Ŝ = M , i.e. M (σw) = Ĉ Âw Bσ for all w ∈ Σ * .

H M α,β , H M σ,α,β , H M α,σ , H M β ((

Correlation

analysis: finding an LPV-SSA representation of y d

In this section, we describe an adaptation of the correlation analysis (CRA) method [START_REF] Cox | Estimation of LPV-SS models with static dependency using correlation analysis[END_REF][START_REF] Cox | Towards efficient maximum likelihood estimation of LPV-SS models[END_REF]

(w) =    1 √ p w E[y(t)(z u w (t)) T ]Λ -1 u ∀w ∈ Σ + E[y(t)u T (t)]Λ -1 u w = (20)
where we recall from Assumption 2, Λ u = var(u).

It

(σ j v j u i ) ki,lj , Ψ u,y (σu i ) ki,l , Ψ u,y (σ j v j ) r,j , i, j = 1, . . . , n x , l = 1, . . . , n u , r = 1, . . . , n y , Ψ u,y ( ).

Covariance realization algorithm

In this section, we adapt the realization algorithm from Mejari and Petreczky (2019a)

P i σ1 ( Âs σ1 ) T + Kσ1 Qi σ1 KT σ1 Qi σ = p σ E[z y s σ (t)(z y s σ (t)) T ] -Ĉs P i σ ( Ĉs ) T Ki σ = Ĝσ √ p σ -Âs σ P i σ ( Ĉs ) T Qi σ -1 (22) with P 0 σ = 0. Moreover, E[e s (t)(e s (t)) T µ 2 σ (t)] = Qσ = lim i→∞ Qi σ , E[x(t)x T (t)µ 2 σ ] = Pσ = lim i→∞ P i σ for all σ ∈ Σ.
The proof of Lemma 5 can be found in [START_REF] Petreczky | Realization theory for a class of stochastic bilinear systems[END_REF]; Mejari and Petreczky (2019a). From Lemma 5, it follows that we can use the basis reduced Kalman-Ho realization algorithm Algorithm 2, as described in Algorithm 3, in order to compute LPV-SSA representation of y s .

Algorithm 3 Realization of y s : Computing an LPV-SSA representation of y s using covariances and n-selection. Input: (n, n y , n y )-selection ᾱ, β of the form of the form (15); Ψ y s (σ j v j u i ) ki,lj ; Ψ y s (σu i ) ki,l ; Ψ y s (σ j v j ) r,j , i, j = 1, . . . , n; l, r = 1, . . . , n y , {E[z y s σ (t)(z y s σ (t)) T ]} σ∈Σ ; number of maximal iterations I > 0. (2) The process (y, u, {µ w } w∈Σ + ) is ergodic and there exist sample paths y : Z → R ny , u : Z → R nu and µ : Z → R nµ of the processes y, u and µ respectively

Construct the matrices H

  requirements for the other pages Paper size this page A4

  ) with the tuple S = ({A σ , B σ } σ∈Σ , C, D). The number n x is called the dimension of S . The sub-Markov parameters of S = ({A σ , B σ } σ∈Σ , C, D) are the values of the map M S : Σ * → R ny×nu , such that for all w ∈ Σ * , M S (w) = CA s B σ , w = σs, σ ∈ Σ, s ∈ Σ * to M S as the sub-Markov function of the deterministic LPV-SSA representation of S . From Petreczky et al. (2017) it then follows that two deterministic LPV-SSA representations S 1 , S 2 have the same input-output behavior, if and only if their sub-Markov parameters are equal, i.e., M S1 = M S2 . Moreover, the sub-Markov parameters can be determined from the input-output behavior. Below we recall from Cox et al. (2018) an adaptation of this Ho-Kalman-like algorithm, which uses sub-Markov parameters to compute a deterministic LPV-SSA representation. In order to present the algorithm, we present the notion of n-selection. Let us define the set Σ n as the set of all words w ∈ Σ * of length less than or equal to n, i.e., Σ n = {w ∈ Σ * | |w| ≤ n}. Definition 5. (Selection). We define (n, n y , n u )-selection as a pair (α, β) such that 1

Algorithm 1

 1 Deterministic realization: Matrix computations using Hankel matrices and n-selection Input: (n, n y , n u )-selection (α, β); Hankel matrix H M α,β (16); shifted Hankel-matrix H M σ,α,β (17); Hankel matrices H M α,σ and H M β defined in (18)-(19) respectively, and M ( ). 1. Compute the matrices Âσ

  M with Ψ u,y in (16)-(19).

  = Ψ u,y ( ). Denote by ({ Âσ , Bσ } σ∈Σ , Ĉ, D) the matrices returned by Algorithm 1. Output: Matrices ({ Âσ , Bσ } σ∈Σ , Ĉ, D) clear from Lemma 4 and Lemma 3 that Algorithm 2 is correct.

  = I ny . Denote by ({ Âσ , Ĝσ } σ∈Σ , Ĉ, I ny ) the matrices returned by Algorithm 1.3. Define Âsσ = 1 √ pσ Âσ , σ ∈ Σ, Ĉs = Ĉ. 4. Compute { Ki σ , Qi σ , P i σ } I i=1 using the recursion (22). Output: Matrices ({ Âs σ , Ĝσ , KI σ , QI σ , P I σ } σ∈Σ , Ĉs )It is clear from Lemma 5 and Lemma 3 that Algorithm 3 is correct.6. IDENTIFICATION ALGORITHMIn this section, we formulate an identification algorithm based on stochastic realization Algorithms 2-3 and selections, for N -length observation sequence of outputs, inputs and scheduling signals, as detailed in Algorithm 4. Intuitively, the main idea behind Algorithm 4 is to estimate the covariances Ψ u,y , Ψ y s and E[z y σ (t)(z y σ (t)) T ] from the observed data and then apply Algorithms 2-3 to the thus estimated covariances. More specifically, the following assumptions are made: Assumption 3. (1) The n x -selection pair (α, β) and ᾱ, β are such that rank H Ψu,y α,β = n x , rank H Ψ y s ᾱ, β = n x , where n x is the state-space dimension of a minimal LPV-SSA realization of y.

  Let us define the map Ψ u,y : Σ * → R ny×nu as follows Ψ u,y

	Margin requirements for the other pages
	Paper size this page A4
	0.986 in
	25 mm
	43 pt	43 pt
	0.597 in	0.597 in
	15.2 mm	15.2 mm
	for finding
	a stationary LPV-SSA representation of y d with noise
	process u.
	71 pt
	0.986 in
	25 mm

  turns out that if y has a stationary LPV-SSA representation with input u, then Ψ u,y is the sub-Markov function of a deterministic LPV-SSA representation. Lemma 4. Assume that ({A σ , B σ }, C, D, x, u) is a stationary LPV-SSA representation (without inputs, Definition 2) of y d . Then Ψ u,y in (20) equals the sub-Markov function M S (14) of the deterministic LPV-SSA representation S = ({A σ , B σ } σ∈Σ , C, D). Conversely, if Ŝ = ({ Âσ , Bσ } σ∈Σ , Ĉ, D) is a deterministic LPV-SSA representation such that its sub-Markov function M Ŝ equals Ψ u,y and it is minimal dimensional among such deterministic LPV-SSA representations, then ({ Âσ , Bσ }, Ĉ, D, x, u) is a stationary LPV-SSA representation (without inputs) of y d .

	The proof of Lemma 4 is presented in Mejari and Petreczky
	(2019b).
	Hence, we can adapt the basis reduced Ho-Kalman
	realization algorithm as described in Algorithm 2. It is
	Algorithm 2 Realization of y

d : Computing an LPV-SSA representation of y d using covariances and n-selection. Input: (n, n y , n u )-selection (α, β) of the form (15); Ψ u,y

  to estimate the stochastic part (7) of a LPV-SSA representation. Define the covariance sequence Ψ y s : Σ * → R ny×ny , where Ψ y s ( ) = I ny , and for all w ∈ Σ + , Lemma 5. If S = ({ Âσ , Ĝσ } σ∈Σ , Ĉ, I ny ) is a minimal dimensional deterministic LPV-SSA representation such that M S = Ψ y s , then ({ Âs σ , Kσ } σ∈Σ , Ĉ, I ny , x, e s ) is a stationary LPV-SSA representation of y s in forward innovation form, where Âs σ =

					1 √ pσ Âσ , Ĉs = Ĉσ ,
		Kσ = lim i→∞	Ki σ , and { Ki σ } σ∈Σ,i∈N satisfies the following
		recursion
		P i+1 σ	=	p σ	Âs σ1
			σ1∈Σ
	Ψ y s (w) = E[y s (t)(z y s w (t)) T ]	(21)
	If y s has a stationary LPV-SSA representation, then
	Ψ y s is a sub-Markov function of a suitable deterministic
	LPV-SSA representation, Petreczky and Vidal (2018);
	Mejari and Petreczky (2019a). Conversely, from a
	deterministic LPV-SSA representation, sub-Markov
	function of which equals Ψ y s a stationary LPV-SSA
	representation can be computed.

Margin requirements for the other pages Paper size this page A4 such that {y(t), u(t), {µ σ (t)} σ∈Σ } N t=1 is observed and the following holds: for all w ∈ Σ * , σ ∈ Σ,

Then for all w ∈ Σ * , σ ∈ Σ, Ψ u,y (w) = lim

where, for all

Lemma 6. (Consistency). With the Assumption 3 the result of Algorithm 4 satisfies the following:

where

We can then view Ψ N y s (w) as an approximation of Ψ y s (w), and

We could modify Algorithm 4 by replacing ( 23) with (25). We conjecture that Lemma 6 will remain true for the modified algorithm.

NUMERICAL EXAMPLE

In this section, we present a numerical example to test the effectiveness of our algorithm. All computations are carried out on an i5 1.8-GHz Intel core processor with 8 GB of RAM running MATLAB R2018a.

Algorithm 4 Identification of stochastic LPV-SSA from observed data. Input: Observations sequence {y(t), u(t), {µ σ (t)} σ∈Σ } N t=1 , and n x -selection (α, β) and ᾱ, β ; {p σ } σ∈Σ , Λ u , maximum number of iterations I > 0.

1. Compute empirical covariances Ψ N u,y (w), for every w ∈ Σ + , such that w = ivu or w = ivσu or w = iv or w = σu for some words v, u ∈ Σ * , σ ∈ Σ, (u, k) ∈ α, (i, v, l) ∈ β for some k = 1, . . . , n y , l = 1, . . . , n u , ∀w ∈ Σ + . 2. Run Algorithm 2 with empirical covariances Ψ N u,y (w), instead of the covariances Ψ u,y . Denote the result returned by Algorithm 2 by S = ({ Ãd σ , Bd σ } σ∈Σ , Cd , Dd ). 3. Compute approximate covariances:

for all σ ∈ Σ and for every w ∈ Σ + , such that w = ivu or w = iv or w = iu or w = ivσu for some words var(yi-ŷi) var(yi) , 0 × 100%, where ŷi denotes the simulated one-step ahead model output and ȳi denotes the sample mean of the output over the validation set.

The LPV-SSA representation in form ( 1) is used for data generation with following matrices:

which corresponds to state-dimension n x = 3, output dimension n y =1, and scheduling signal dimension n µ =2 with Σ = {1, 2}.

Training and noise free validation output sequences of length N = 100000 and N val = 100000, respectively, are generated using a white-noise input process u with uniform distribution U(-1.5, 1.5) and an independent scheduling signal process µ = [µ 1 µ 2 ] such that µ 1 (t) = 1 and µ 2 (t) is a white-noise process with uniform distribution U(-1.5, 1.5). This corresponds to the parameter values {p σ } σ∈{1,2} to be p 1 =E[µ 2 1 (t)]=1 and p 2 =E[µ 2 2 (t)]=0.75. The standard deviation of the white Gaussian noise e corrupting the training output is 1, i.e., e ∼ N (0, 1). This corresponds to the Signal-to-Noise Ratio SNR = 10 log We run the version of Algorithm 4 explained in Remark 1, with I = 50 iterations and with the following n-selection pairs (α, β) and (ᾱ, β), with n = 3, α = {( , 1), (1, 1), (21, 1)}, β = {(2, , 1), (1, 2, 1), (2, 21, 1)}, ᾱ = {( , 1), (1, 1), (21, 1)}, β = {(1, , 1), (1, 2, 1), (1, 21, 1)}, which are used to choose corresponding entries of the Hankel matrices. The mean time taken to run the algorithm is 1.55 sec. The validation result using one-step ahead predicted outputs ŷ are reported in Table 1, and true vs estimated sub-Markov parameters are reported in Table 2. The results show a good match between estimated model output w.r.t. true system output.

CONCLUSION

In this paper, we formulated a realization algorithm and an efficient identification algorithm for stochastic LPV-SSA representations with inputs, by combining correlation analysis method with a stochastic realization based identification algorithm. The proposed algorithm provides a computationally efficient alternative to the parametric subspace approaches avoiding the curse of dimensionality.