Identification of Multivariable Canonical
State-Space Representations
L. Bako ∗ G. Merce`re ∗∗ and S. Lecœuche ∗
∗ Ecole des Mines de Douai, De´partement Informatique et
Automatique, 941, rue Charles Bourseul, 59508 Douai Cedex, France
(e-mail: bako@ensm-douai.fr, lecoeuche@ensm-douai.fr)
∗∗ Universite´ de Poitiers, Laboratoire d’Automatique et d’Informatique
Industrielle, 40 avenue du recteur Pineau, 86022 Poitiers, France
(e-mail: guillaume.mercere@univ-poitiers.fr)
Abstract: The paper introduces a method to estimate a linear state-space model from its
input-output data. Its starting point is the conversion of the state-space model into an input-
output ARMAX-like model. In order to conveniently estimate the parameters of this ARMAX
model, it is necessary, given the structural characteristics of the initial state-space model, to
assume that the dynamic state matrix A is nonderogatory. Under this condition, it is shown
that the input-output model and the state-space model from which it is derived have the same
order. Based on this important result in linear system realization theory, a canonical state-
space representation is finally obtained through a simple construction technique. The developed
method can be regarded as an alternative to the standard subspace methods particularly in
situations where the state coordinates basis need to be maintained constant with respect to the
input-output realization used for the estimation.
Keywords: Least-squares algorithm ; State-space representation ; Canonical form ; ARMAX
model ; Multi-input/multi-output systems.
1. INTRODUCTION
During the last two decades, interesting developments
have been achieved in linear system identification thanks
the so-called subspace-based methods (Katayama, 2005).
One of the main advantages of this approach is that it
gives directly access to minimal state-space representation
without requiring the use of canonical parameterizations
and the determination of structural indices demanded by
classic prediction error methods when multivariable sys-
tems are studied. Although the subspace-based methods
have proved their efficiency from a practical point of view
(Mevel et al., 1999), (Favoreel et al., 2000), (Lovera, 2003),
most of the techniques developed until now suffer from
two main drawbacks which can be crucial in particular
cases, e.g, when hybrid systems or on-line identification
framework are involved. More precisely,
• the user of the common subspace algorithms has
difficulties to fix the basis in which the state-space
models are estimated. In fact, this choice is realized
during the reduction step (when the system order
is estimated) by using generally a Singular Value
Decomposition (SVD) which gives rise to particular
fully-parametrizations. But the reproducibility of the
basis is not ensured when such a tool is employed.
This characteristic can be problematic when hybrid
systems such as switched or LPV processes are identi-
fied since the different sub-models must be expressed
in coherent bases. Otherwise, the input-output behav-
ior of the global model may be distorted (To´th et al.,
2007).
• the SVD, interesting for its computational robustness,
is numerically cumbersome. Thus, this tool is diffi-
cult to use in a recursive estimation framework and
alternatives must be developed (Lovera et al., 2000),
(Merce`re et al., 2008) to overcome this problem when
it is necessary to update the model in real time.
To sort out these two problems, a new SVD-free algorithm
for the identification of linear MIMO state-space models is
suggested in this paper. It is based on an initial estimation
step of an ARMAX model using traditional least squares
algorithms followed by an algebraic reconstruction of a
canonical state-space representation from the parameters
calculated during the first phase. This approach is inspired
by the works of Stoica and Jansson (2000) where it
is shown that canonical transfer function approach can
be seen as an interesting alternative to the subspace
methods even when MIMO systems are considered. The
introduction of least squares optimization techniques is an
important asset of the developed method since it avoids
the use of SVD and make its adaptation easier for on-line
scenarii. Furthermore, as shown afterwards, the estimated
state-space matrices will be given in a fixed coordinate
basis. Thus, this algorithm can be easily extended when
composite systems are involved.
The outline of the paper is as follows. Section 2 is con-
cerned with the presentation of the model, the relative
assumptions and a statement of the identification problem.
In Section 3, the state-space model is transformed into

an ARMAX-like model whose parameters are estimated
under an appropriate assumption on the A-matrix. Given
this latter model, an original but quite simple technique
to construct a state-space model is suggested in Section 4.
Conditions that guarantee equality of orders between the
initial state-space model and its corresponding ARMAX-
like model are also derived.
2. PROBLEM FORMULATION AND NOTATIONS
Consider the following finite dimensional, causal and linear
time-invariant state-space model
x(t + 1) = Ax(t) +Bu(t) +w(t) (1a)
y(t) = Cx(t) +Du(t) + v(t) (1b)
where u(t) ∈ Rnu , y(t) ∈ Rny and x(t) ∈ Rn are re-
spectively the input, the output and the state vectors.
The process noise sequence is w(t) ∈ Rn and the output
measurement disturbances are collected in v(t) ∈ Rny .
(A,B,C,D) are the system matrices relatively to a certain
coordinate basis of the state-space. It will be assumed
that {x(t)}, {u(t)}, {y(t)}, {w(t)} and {v(t)}, t ∈ Z, are
all ergodic and (weakly) stationary stochastic processes
(Ljung, 1999), the system is asymptotically stable, control-
lable and observable and there is no feedback from y to u.
The considered identification problem can then be stated
as follows: given realizations {u(t)}Nt=1 and {y(t)}
N
t=1 of
the input and output processes generated by a system of
the form (1) on a finite but sufficiently wide time horizon
N , estimate the matrices (A,B,C,D) up to a similarity
transformation.
To reach this goal, two steps will be considered. Firstly,
an input-output ARMAX model will be identified from
the input-output measurements by using a classic least
squares or instrumental variable method. Then, under the
assumption that the matrix A is nonderogatory (Horn and
Johnson, 1990), a canonical state-space form will be recon-
structed from the estimated ARMAX model parameters.
3. ARMAX MODEL ESTIMATION
Consider the output equation (1b). A direct way to get an
input-output representation of the system from this state-
space form is to get rid of the state. For that, let us start
by expressing the state as a linear combination of the past
inputs, outputs and a prior value of the state. Particularly,
for t > n, it is easy to show that
x(t) = Anx(t − n) +∆nun(t− n) + Anwn(t − n) (2)
with 1
un(t − n) =
[
u⊤(t − n) · · · u⊤(t − 1)
]⊤ ∈ Rnnu
∆n =
[
An−1B · · · AB B
]
∈ Rn×nnu
An =
[
An−1 · · · A In
]
∈ Rn×n2 .
Substituting Eq. (2) into Eq. (1b), we get
y(t) = CAnx(t − n) +C∆nun(t − n)
+CAnwn(t− n) +Du(t) + v(t). (3)
Now, by using the Cayley-Hamilton theorem (Horn and
Johnson, 1990), we know that An = −a0In − a1A −
1 yn(t − n) and wn(t− n) are built in a similar way.
· · · − an−1An−1 where the {ai}i∈[0,n−1] elements are the
coefficients of the characteristic polynomial of A, i.e.,
pA(z) = a0 + a1z + · · ·+ an−1zn−1 + zn.
Thus,
CAnx(t − n) = −
(
a⊤ ⊗ Iny
)
Γnx(t − n) (4)
where
a⊤ = [a0 a1 · · · an−1] (5)
Γn =
[
C⊤ (AC)⊤ · · ·
(
CAn−1
)⊤
]⊤
(6)
and with ⊗ the Kronecker product (Horn and Johnson,
1990). By using now the recurrent relations (1), it is
straightforward to verify that the stacked vectors yn(t−n),
un(t − n), vn(t − n) and wn(t − n) are related as follows
yn(t − n) = Γnx(t − n) +Hnun(t − n)
+Gnwn(t − n) + vn(t − n) (7)
with
Hn =



D 0 ··· 0
CB
. . . ...
... . . . . .. 0
CAn−2B ··· CB D


 , Gn =



0 0 ··· 0
C
. .. ...
... . .. . .. 0
CAn−2 ··· C 0


 .
Eq. (7) allows us to express Γnx(t − n) as a combination
of past inputs, outputs and disturbances, i.e.,
Γnx(t − n) = yn(t − n)−Hnun(t − n)
−Gnwn(t − n)− vn(t − n). (8)
Finally, gathering equations (3), (4) and (8), the state
vector is eliminated in the output equation and we get
y(t) = −
(
a⊤ ⊗ Iny
)
yn(t− n) +Fun+1(t− n) + η(t) (9)
with
F =
[(
a⊤ ⊗ Iny
)
Hn +C∆n D
]
∈ Rny×(n+1)nu
η(t) =
((
a⊤ ⊗ Iny
)
Gn +CAn
)
wn(t− n)
+
(
a⊤ ⊗ Iny
)
vn(t − n) + v(t).
Note that F contain the coefficients of the polynomial ma-
trix 2 that defines the numerator of the transfer function of
(1) and a is the coefficient vector of the denominator. The
model (9) is in fact quite more complex than a standard
ARMAX since it contains at the same time a measurement
and a process noise. It is also interesting to note that the
noise part of this ARMAX model is composed of particular
combinations of the AR and X deterministic parameters.
Thus, as soon as a and F are consistently estimated, a
model of the stochastic part can be built.
In order to simplify the estimation of a and F, it is
interesting to notice that Eq. (9) can be rewritten as a
classic linear regression. For that, introduce the identity
vec (AXB) =
(
B⊤ ⊗A
)
vec (X) where vec (•) is the
vectorization operator (Horn and Johnson, 1990). Then,
y(t) =
[
−Φyn(t − n) u⊤n+1(t − n)⊗ Iny
]
[
a
vec (F)
]
+ η(t)
(10)
with Φyn(t − n) = [y(t − n) · · · y(t − 1)] ∈ Rny×n. Thus,
a and F can be estimated using a least squares algorithm
such as the instrumental variable method (Ljung, 1999)
when the data matrix
2 For an explicit expression of this polynomial matrix, we refer to
the proof of Theorem 1 given in Appendix B.

Φ(n) =



Φyn(1) u⊤n+1(1)⊗ Iny
...
...
Φyn(N − n) u⊤n+1(N − n)⊗ Iny



has full column rank.
4. FROM ARMAX MODEL TO MINIMAL
CANONICAL STATE-SPACE REPRESENTATION
4.1 Condition to get ARMAX and state-space models of
the same order
The problem of parameterizing state-space models and the
equivalence between state-space representation and MIMO
canonical transfer functions is not new (see, e.g. (Kailath,
1980; McKelvey, 1995)). One of the main problems is
the relation between the order nIO of the I/O model
and the one of the corresponding minimal state-space
representation. A glance at Eq. (9) and the way it is built
from the state-space model (1) may tend to suggest that
both representations have the same order. Unfortunately,
this is not the case for all the MIMO systems (see (Stoica
and Jansson, 2000) for a counterexample). Indeed, we have
in general,
nIO ≤ nSS ≤ nIO · ny,
where nSS = n is the order of the state-space model (1).
To see this fact more clearly, let us think of the order as
the smallest integer m such that y(t) can be expressed in
function of ym(t−m) and um+1(t−m) as in Eq. (9). Then,
intuitively, the order nIO of the input-output model (9)
should be equal to the degree ν of the minimal polynomial
of A. This remark follows essentially from Eq. (4) which
suggests that nIO should be the smallest integer such that
AnIO can be expressed in function ofAj , j = 0, . . . , nIO−1,
where A is taken from a minimal realization of (1).
On the other hand, if we consider the input-output transfer
matrix of (1) to be TF(z) = L (z)/pA(z), with
L (z) = Cadj(zI −A)B+ pA(z)D
=



f11(z) · · · f1nu(z)
...
...
fny1(z) · · · fnynu(z)


 ,
then the order nIO of Eq. (9) corresponds to
n− deg
{
gcd
(
{fij(z)}j∈[1,nu]i∈[1,ny ] , pA(z)
)}
.
In fact, minimality of (A,B,C,D) is not enough to
guarantee coprimeness between the polynomials fij(z) and
pA(z). This is formally stated by the following theorem
which relates the order of model (9) to that of (1).
Theorem 1. Assume that the model (1) is minimal, that is,
observable and reachable. Then the order nIO of the input-
output model (9) is equal to the degree ν of the minimal
polynomial of A, i.e., nIO = deg(mA(z)) = ν ≤ n.
For the proof, see Appendix A.
We know from Theorem 1 that, in general, the orders of
the models (1) and (9) do not coincide, even when the
model (1) is minimal. Therefore, even though it is assumed
that the order n is known a priori, it is compulsory to
determine nIO before being able to estimate a and F from
(10). Indeed, consistent estimation of these parameters
requires that the matrix Φ(n) has full column rank. No
matter how rich is the excitation input, it can be easily
shown that rank(Φ(n)) ≤ nIO + (n+ 1)nuny in the noise-
free case (i.e., when η(t) = 0 in Eq. (9)). It follows that
it is necessary to satisfy nIO = n in order to have a full
column rank matrix Φ(n).
To get equality between both orders, a special property on
the matrix A in Eq. (1) is required as formulated by the
following corollary of Theorem 1.
Corollary 2. Assume that the model (1) is minimal. Then
nIO = nSS = n if and only if the matrix A is nonderoga-
tory. 3
In order to explicitly reconstruct a state-space representa-
tion from the estimated ARMAX parameters, it is essential
to find conditions ensuring an equality relation between
nIO and nSS.
4.2 Reconstruction of a minimal canonical state-space
model
Once the parameters of the model (9) are estimated, the
reconstruction of a canonical state-space model can be
considered. Before, the following lemma must be stated.
Lemma 3. If the system (1) is observable, then there exists
a matrix Ψ of the form
Ψ =






γ⊤ 0 · · · 0
0 . . .
...
... . . . 0
0 · · · 0 γ⊤






∈ Rn×nny , (11)
with γ ∈ Rny , and verifying
rank (ΨΓn) = n (12)
where Γn is the observability matrix of order n, if and only
if the matrix A is nonderogatory.
A proof of this lemma can be found in (Bako et al.,
2009). Note that the condition (12) in Lemma 3 can be
interpreted as follows. There exists a γ such that if we let
ya(t) = γ⊤y(t) be a linear combination of all the outputs,
then the state x(t) of system (1) is entirely observable
from ya(t). This further means that the system (1) and
the MISO system with input u(t) and output ya(t) ∈ R
have the same state. Lemma 3 tells us that this is in fact
possible if and only if the matrix A is nonderogatory 4 . In
view of that lemma, if, e.g, γ ∈ Rny is randomly generated
from a uniform distribution, then the required condition
(12) is fulfilled with probability one (Bako et al., 2009).
By assuming now that the matrix A is nonderogatory
(which ensures the equality nIO = nSS), a minimal state-
space realization can be readily obtained as explained in
the following proposition:
Proposition 4. Assume that the system (1) is minimal
and that the matrix A is nonderogatory. Consider a
matrix Ψ defined as in (11) and satisfying condition
3 A square matrix is said to be nonderogatory if its minimal
polynomial coincides with its characteristic polynomial (Horn and
Johnson, 1990).
4 In fact, for any SISO or MISO system to be observable the A-
matrix needs to be nonderogatory.

(12). Then a minimal canonical state-space representation
can be obtained from the estimated parameters a⊤ =
[a0 · · · an−1] and F of the ARMAX model (9) as:
A =




0 1 · · · 0
...
... . . . 0
0 0 · · · 1
−a0 −a1 · · · −an−1




, (13)
D = F (:, nnu + 1 : (n + 1)nu) , (14)
B = ΨQ, (15)
C = rowny (Q)Ω⊤n
(
ΩnΩ⊤n
)−1 , (16)
where
Q =
(
K⊗ Iny
)−1 colnu
(
F(:, 1 : nnu)− a⊤ ⊗D
)
,
Ωn =
[
B AB · · · An−1B
]
∈ Rn×nnu ,
K =




a1 · · · an−1 1
... 1 0
an−1
1 0 · · · 0




∈ Rn×n,
rowny
([ X1
...
Xn
])
= [X1 ··· Xn ] , colnu ([Y1 ··· Yn ]) =
[ Y1
...
Yn
]
,
with Xi ∈ Rnyו and Yi ∈ R•×nu , i ∈ [1, n].
For the proof, see appendix B.
As opposed to classic techniques leading to canonical forms
for multivariable systems (Luenberger, 1967), (Mayne,
1972), (Kailath, 1980), the developed approach does not
need any observability or controllability indices or an ini-
tial step leading to a non minimal realization estimate
followed by a reduction procedure. Here is proposed a
relatively straightforward method which gives access to
minimal MIMO state-space representation from a stan-
dard vector difference equation.
5. SIMULATION EXAMPLE
In order to illustrate the performances of the algorithm,
the following system is considered
A =
[
−0.3814 0.6134 −0.3495
0.4044 −0.0624 −0.7160
−0.5787 −0.5476 −0.1790
]
, B =
[
0.8736 0
0 −0.3881
0 0
]
,
C =
[0.9397 0 1.1787
0 0 −1.3274
]
, D =
[0.5463 −0.5293
0 −2.4003
]
.
The input is a zero-mean white Gaussian noise with unit
variance. The disturbances v and w are both zero-mean
white Gaussian noise set such that the noise signal ratio
equals 20dB in both cases. A Monte Carlo simulation with
100 realizations of noise and input is carried out.
This new algorithm is compared with three other methods:
• another SVD-free identification method developed in
(Bako et al., 2008) where the subspace identification
problem is recasted into a simple least squares prob-
lem by introducing a particular user-defined matrix
which does not change the rank of the extended ob-
servability matrix when multiplying this latter matrix
on the left hand side,
• the traditional PO MOESP algorithm (Verhaegen,
1994),
• an off-line version of the propagator method sug-
gested initially in (Merce`re et al., 2003) in a recur-
sive form (see also (Merce`re et al., 2008)) where a
particular linear operator is introduced to rewrite
the subspace-based identification problem into a
quadratic optimization problem without constraint.
For all these methods, the instrumental variable used to
reduce the effect of the noise is built from the inputs of the
system and the outputs of a simulated biased auxiliary
model obtained with the ordinary MOESP algorithm
(Verhaegen, 1994) applied on noisy data.
−1 −0.5 0 0.5 1
−0.2
−0.1
0
0.1
0.2
0.3
Real part
Im
ag
in
ar
y
pa
rt
−1 −0.5 0 0.5 1
−0.2
−0.1
0
0.1
0.2
0.3
Real part
Im
ag
in
ar
y
pa
rt
−1 −0.5 0 0.5 1
−0.2
−0.1
0
0.1
0.2
0.3
Real part
Im
ag
in
ar
y
pa
rt
−1 −0.5 0 0.5 1
−0.2
−0.1
0
0.1
0.2
0.3
Real part
Im
ag
in
ar
y
pa
rt
Method in (Bako et al., 2008)
Method in (Merce`re et al., 2008)PO-MOESP
Our algorithm
Fig. 1. Estimated (o) and real poles (+) of the system.
Figure 1 compares the position of the estimated poles of
the models given by the four aforementioned algorithms
with respect to the real system. A first conclusion is that
all these methods are quite accurate since the estimated
poles are centered on the true ones and the variances of
the estimates are relatively small. We can however notice
that, on this example, the new algorithm leads, without
requiring the use of an SVD, to similar results as the PO
MOESP algorithm. Therefore, the approach developed in
this article seems to be an interesting alternative to the
classic subspace-based techniques.
6. CONCLUSION
In this paper, a simple SVD-free method has been pro-
posed to get a state-space model, given a collection of its
input-output measurements. In comparison with existing
subspace identification techniques, the main motivation of
this new method is to prevent the state representation
basis from shifting according to the estimation data. This
feature is highly desirable in, e.g., recursive identification
of state-space models or in the estimation of composite
models with linear submodels. From a realization theory
point of view, this work highlights the behavior of certain
types of linear MIMO systems that can indeed be handled
as MISO systems in terms of state dynamics. Future re-
search will focus on generalizing this procedure to more
arbitrary models than the particular cases where the A-
matrix is nonderogatory.

Appendix A. PROOF OF THEOREM 1
We first state the following lemma.
Lemma 5. Let A ∈ Rn×n and B = [b1 · · · bm] ∈ Rn×m
be some arbitrary matrices. Let mA(z) be the mini-
mal polynomial of A and define
[
An−1B · · · AB B
]
=
∆n(A,B) ∈ Rn×nm,
A¯ =



A
. . .
A


 ∈ Rnm×nm, b¯ =



b1
...
bm


 ∈ Rnm.
Then the following holds. If rank
(
∆n(A,B)
)
= n, then
rank
(
∆n(A¯, b¯)
)
= ν, where ν = deg(mA(z)).
Proof. Assume that rank
(
∆n(A,B)
)
= n. We then need
to show that the rank of the matrix ∆n(A¯, b¯) is equal to
the degree ν of the minimal polynomial of A. We proceed
by contradiction. Assume that there exists j < ν such
that A¯j b¯ is a linear combination of the vectors A¯ib¯ ,
i = 0, . . . , j − 1. Then, one can find a set of scalars
α0, · · · , αj−1 not all zero satisfying A¯j b¯ = αj−1A¯j−1b¯ +
· · · + α0b¯. By Denoting with p the polynomial of degree
j defined by p(z) = zj − αj−1zj−1 − · · · − α0, we get
p(A¯)b¯ = 0. This is equivalent to p(A)bi = 0, i = 1, . . . ,m.
As a consequence, p(A)B = 0 and p(A)∆n(A,B) = 0
because, as a polynomial of A, p(A) commutes with
any power of A. Since ∆n(A,B) has full row rank, it
follows that p(A) = 0. Therefore, if we call mA(z) the
minimal polynomial of A, then, by definition of mA(z)
we get that mA(z) divides p(z). This contradicts the fact
that deg {p(z)} = j < deg {mA(z)} = ν. In conclusion,
rank
(
∆n(A¯, b¯)
)
= ν.
From this lemma, it can be concluded that under control-
lability assumption of system (1), rank
(
∆n(A¯, b¯)
)
= n if
and only if the A-matrix is nonderogatory. For the proof of
Theorem 1, we will need also the following Theorem given
in (Barnett, 1973):
Theorem 6. Let
bi(z) = bi1zn−1 + · · · bin, i = 1, · · · , p,
pA(z) = a0 + a1z + · · ·+ an−1zn−1 + zn.
Then the degree of the gcd of pA(z), b1(z), · · · , bp(z) is
equal to n− rank (Γn (Ac,M)) with
M =



b1n · · · b11
...
...
bpn · · · bp1


 , Ac =




0 1 · · · 0
...
... . . .
0 0 · · · 1
−a0 −a1 · · · −an−1




.
With this material, Theorem 1 can be proved as follows:
Proof. The input-output transfert matrix is given by
C(zI −A)−1B+D = Cadj(zI −A)B+ pA(z)DpA(z)
,
where adj(·) is the cofactor matrix of (·) and pA(z) =
det(zI −A) is the characteristic polynomial of A. Let us
introduce the following notation
Cadj(zI−A)B = [lij(z)]j∈[1,nu]i∈[1,ny] ,
where lij(z), the numerator polynomial that corresponds
to the transfert between the i-th output and the j-th input.
To prove Theorem 1, we need to show that the degree of
gcd
(
{lij(z)}j∈[1,nu]i∈[1,ny ] , pA(z)
)
is equal to n − ν, with ν = deg(mA(z)). To proceed, we
start by pointing out that, as far as we are concerned with
seeking the common factors between the polynomials lij(z)
and pA(z), we can, without loss of generality, assume that
D = 0.
Now we compute the coefficients of the polynomials lij(z).
The expressions of these coefficients can be directly ob-
tained from the term F of Eq. (9) (i.e., the matrix that
multiplies un(t−n)). At this step, it is important to notice
that
F = (a⊤ ⊗ Iny )Hn +C∆n = C∆n(L⊗ Inu),
with
L =





1 0 · · · 0
an−1 1
...
... . . . . . . 0
a1 · · · an−1 1





∈ Rn×n.
Let bj be j-th column of B and c⊤i be the i-th row of C.
The n coefficients of the n−1 degree polynomial lij(z) are
given by the following row vector
lij = c⊤i
[
An−1bj · · · Abj bj
]
︸ ︷︷ ︸
∆n(A,bj)
L.
Let
Nj =



l1j
...
lnyj


 = C∆n(A,bj)L
and then
N =



N1
...
Nnu


 =



C∆n(A,b1)
...
C∆n(A,bnu)


L ∈ Rnny×n
as a concatenation of all the coefficient vectors of the
polynomial matrix Cadj(zI−A)B. According to Theorem
6, we just need to show that
rank (P) = ν,
where
P =




N
NAc
...
NAn−1c




and Ac a companion matrix (see Eq. (13)). Notice that
the matrix P is equal to
P1 =












C∆n(A,b1)L
C∆n(A,b1)LAc
...
C∆n(A,b1)LAn−1c...
C∆n(A,bnu)L
C∆n(A,bnu )LAc
...
C∆n(A,bnu)LAn−1c












=













C∆n(A,b1)
C∆n(A,b1)AL
...
C∆n(A,b1)An−1L...
C∆n(A,bnu)
C∆n(A,bnu)AL
...
C∆n(A,bnu)An−1L













L,
up to some rows permutations, where
AL = LAcL−1 =




−an−1 1 · · · 0
...
... . . .
...
−a1 0 · · · 1
−a0 0 · · · 0




.

For any k ≥ 0 and any bj , we have ∆n(A,bj)AkL =
Ak∆n(A,bj). By making use of this fact, it follows that
P1 =



Γn
. . .
Γn


∆n(A¯, b¯)L,
with Γn = Γn(A,C). Since the system is observable, the
matrix that contains Γn in the right hand side of the pre-
vious equation, has full column rank. Therefore, Lemma
5 allows us to conclude that rank
(
P
)
= rank
(
P1
)
=
rank
(
∆n(A¯, b¯)
)
= ν since L is nonsingular.
Appendix B. PROOF OF PROPOSITION 4
The reconstruction of a canonical form from the parame-
ters a and F is in fact only based on shrewd calculations.
For D, it is easy to see that
D = F (:, nnu + 1 : (n + 1)nu) .
Knowing D, it is interesting to notice that
F (:, 1 : nnu)− a⊤ ⊗D
=
(
a⊤ ⊗ Iny
)
Hn +C∆n − a⊤ ⊗D
=
[(
k⊤1 ⊗ Iny
)
ΓnB · · ·
(
k⊤n ⊗ Iny
)
ΓnB
]
where k⊤i is the ith row of K. Thus, by using the operator
colnu defined in Proposition 4 and noticing that K ⊗ Iny
is invertible, we get
ΓnB =
(
K⊗ Iny
)−1 colnu
(
F (:, 1 : nnu)− a⊤ ⊗D
)
.
Introducing matrix Ψ, it holds that
ΨΓnB = TB = B¯
= Ψ
(
K⊗ Iny
)−1 colnu
(
F (:, 1 : nnu)− a⊤ ⊗D
)
since ΨΓn = T ∈ Rn×n can be viewed as a similarity
transformation if rank (ΨΓn) = n. The next step consists
in reconstructing the state matrix A. This can be done by
calculating A¯ = TAT−1 knowing T, i.e.
A¯ = ΨΓnA (ΨΓn)−1
= Ψ



CA
...
CAn


 (ΨΓn)−1 =



γ⊤CA
...
γ⊤CAn


 (ΨΓn)−1 .
Now, using again the Cayley-Hamilton theorem, we have



γ⊤CA
...
γ⊤CAn


 =




0 1 · · · 0
...
... . . .
0 0 · · · 1
−a0 −a1 · · · −an−1







γ⊤C
...
γ⊤CAn−1



= AcΨΓn.
where Ac is a companion matrix (see Eq. (13)). Hence,
A¯ = AcΨΓn (ΨΓn)−1 = Ac.
Finally, knowing A¯, B¯ and ΓnB, it is easy to see that
rowny (ΓnB) = rowny
(
Γ¯nB¯
)
= C¯
[
B¯ · · · A¯n−1B¯
]
which leads to the extraction of C¯.
REFERENCES
L. Bako, G. Merce`re, and S. Lecœuche. A least squares
approach to the subspace identification problem. In
In the 47th IEEE Conference on Decision and Control,
Cancun, Mexico, December 2008.
L. Bako, G. Merce`re, and S. Lecœuche. On-line structured
subspace identification with application to switched
linear systems. International Journal of Control, 2009.
Accepted for publication.
S. Barnett. Matrices, polynomials and linear time invari-
ant systems. IEEE Transactions of Automatic Control,
18:1–10, 1973.
W. Favoreel, B. De Moor, and P. Van Overschee. Subspace
state space system identification for industrial processes.
Journal of Process Control, 10:149–155, 2000.
R. Horn and C. Johnson. Matrix analysis. Cambridge
University Press, 1990.
T. Kailath. Linear Systems. Prentice Hall, Engelwood
Cliffs, 1980.
T. Katayama. Subspace methods for system identification.
Springer, 2005.
L. Ljung. System identification. Theory for the user.
Prentice Hall Information and System Sciences Series,
Upper Saddle River, 2nd edition, 1999.
M. Lovera. Identification of MIMO state space models for
helicopters dynamics. In Proceedings of the 13th IFAC
Symposium on System Identification, Rotterdam, The
Netherlands, August 2003.
M. Lovera, T. Gustafsson, and M. Verhaegen. Recursive
subspace identification of linear and non linear Wiener
state space models. Automatica, 36:1639–1650, 2000.
D. Luenberger. Canonical forms for linear multivariable
systems. IEEE Transactions on Automatic Control, 12:
290–293, 1967.
D. Mayne. A canonical model for identification of mul-
tivariable linear systems. IEEE Transactions on Auto-
matic Control, 17:728–729, 1972.
T. McKelvey. Identification of state space models from
time and frequency data. PhD thesis, Linko¨ping Univer-
sity, Linko¨ping, Sweden, 1995.
G. Merce`re, S. Lecœuche, and C. Vasseur. A new recursive
method for subspace identification of noisy systems:
EIVPM. In Proceedings of the 13th IFAC Symposium
on System Identification, Rotterdam, The Netherlands,
August 2003.
G. Merce`re, L. Bako, and S. Lecœuche. Propagator-based
methods for recursive subspace model identification.
Signal Processing, 88:468–491, 2008.
L. Mevel, L. Hermans, and H. Van der Auweraer. On the
application of subspace-based fault detection methods
to industrial structures. Mechanical Systems and Signal
Processing, 13:823–838, 1999.
P. Stoica and M. Jansson. MIMO system identification:
state-space and subspace approximations versus transfer
function and instrumental variables. IEEE Transactions
on Signal Processing, 48:3087–3099, 2000.
R. To´th, F. Felici, P. Heuberger, and P. Van den Hof.
Discrete-time LPV I/O and state-space representations.
differences of behavior and pitfalls of interpolation. In
Proceeding of the European Control Conference, Kos,
Greece, July 2007.
M. Verhaegen. Identification of the deterministic part of
MIMO state space models given in innovations form
from input output data. Automatica, 30:61–74, 1994.