arXiv:1801.08409v1 [math.NA] 14 Jan 2018

Computation of the State Bias and Initial States
for Stochastic State Space Systems in the General 2-D
Roesser Model Form
José A. Ramos
Nova Southeastern University
College of Engineering and Computing
Department of Engineering and Technology
3301 College Avenue
Fort Lauderdale, FL 33314
Email: jr1284@nova.edu
and
Guillaume Mercère
Université de Poitiers
Laboratoire d’Informatique et d’Automatique pour les Systèmes
2 rue Pierre Brousse, bâtiment B25, TSA 41105
86073 Poitiers cedex 9, France
Email: guillaume.mercere@univ-poitiers.fr
January 26, 2018

Abstract
Recently Ramos and Mercère (2017a) presented a subspace system identification algorithm
for 2-D purely stochastic state space models in the general Roesser form. However, since
cvh ,
the exact problem requires an oblique projection of Y hf projected onto W hp along X
f
"
#
vh
c
cvh , X
cvh } are unknown. In the
where W hp = X ph , this presents a problem since {X
p
f
Yp
above mentioned paper, the authors found that by doing an orthogonal projection Y hf /Y hp ,
ch with a small bias due to the initial
one can identify the future horizontal state matrix X
f

cvh , X
cvh }. Nevertheless, the results on modeling 2-D images
conditions that depend on {X
p
f
cvh , X
cvh }. In this note we delve into the bias
were very good despite lack of knowledge of {X
p
f
term and prove that it is insignificant, provided i is chosen large enough and the vertical and
horizontal states are uncorrelated. That is, the cross covariance of the state estimates xhr,s and
xvr,s is zero, or Phv = 0nx ×nx and Pvh = 0nx ×nx . Our simulations use i = 30. We also present
a second iteration to improve the state estimates by including the vertical states computed
from a vertical data processing step, i.e., by doing an orthogonal projection Y vf /Y vp . In
this revised algorithm we include a step to compute the initial states. This new portion, in
addition to the algorithm presented in Ramos and Mercère (2017a), forms a complete 2-D
stochastic subspace system identification algorithm.

1

Problem Formulation

The general 2-D stochastic Roesser model has the state-space form
h
xhr+1,s = A1 xhr,s + A2 xvr,s + wr,s

(1a)

xvr,s+1

(1b)

=

yr,s =

A3 xhr,s
C1 xhr,s

v
+ A4 xvr,s + wr,s
+ C2 xvr,s + vr,s ,

(1c)

where xhr,s ∈ Rnh , xvr,s ∈ Rnv , and yr,s ∈ Rny denote, respectively, the local horizontal
state, local vertical state, and output vectors at the (r, s)th location of a finite domain
D = {(r, s) | 0 ≤ r ≤ N and 0 ≤ s ≤ M}. The system matrices {A, C}, given by




A1 A2
(2a)
, C = C1 C2 ,
A =
A3 A4

have partitioned dimensions A1 ∈ Rnh ×nh , A2 ∈ Rnh ×nv , A3 ∈ Rnv ×nh , A4 ∈ Rnv ×nv , C1 ∈
h
v
Rny ×nh , and C2 ∈ Rny ×nv . The noise vectors wr,s
∈ Rnh , wr,s
∈ Rnv , and vr,s ∈ Rny are
assumed to be white Gaussian noise processes with mean and joint covariance matrix given,
respectively, by
 h 


0
w


n
×1
r,s
h




 v 
IE  wr,s
=
(3a)
0
 nv ×1 





vr,s
0ny ×1





h
Q
Q
S
w


hh
hv
h
r,s

h
⊤ ⊤ i
⊤


 v 
=  Qvh Qvv Sv  · δr−r′ · δs−s′
IE  wr,s
vr′ ,s′
wrv′ ,s′
wrh′ ,s′





vr,s
Sh⊤ Sv⊤ R
"
#
Q S
=
(3b)
· δr−r′ · δs−s′ ,
S⊤ R

where Qhh ∈ Rnh ×nh , Qhv ∈ Rnh ×nv , Qvh ∈ Rnv ×nh , Qvv ∈ Rnv ×nv , Sh ∈ Rnh ×ny , Sv ∈ Rnv ×ny ,
and R ∈ Rny ×ny , nx = nh + nv is the dimension of the combined system, IE is the expectation
operator, M ⊤ denotes the transpose of M, δk−k′ is the Kronecker delta function, 0m×n denotes
an (m × n) matrix with all its elements equal to zero, and {Q, R, S} are the covariance and
cross-covariance matrices of the noise terms.
The noise and state vectors are uncorrelated with each other, i.e.,
n
h
⊤
⊤ ⊤ io
h
h
v
IE xr,s
= 0nh×(nx +ny ) , ∀ r ′ ≥ r and s′ ≥ s
(4a)
wr′ ,s′
wr′ ,s′
vr′ ,s′
io
h
n
⊤ ⊤
⊤
= 0nv ×(nx +ny ) , ∀ r ′ ≥ r and s′ ≥ s.
(4b)
IE xvr,s
wrv′ ,s′
vr′ ,s′
wrh′ ,s′
Furthermore, the states xhr,s and xvr,s evolve with the following statistical properties: zero
mean
#
#)
"
("
0nh×1
xhr,s
, r = 0, 1, . . . , N and s = 0, 1, . . . , M
(5)
=
IE
xvr,s
0nv ×1
1

and positive definite state covariance matrix
#
#
("
)
"
i
xhr,s h h ⊤
0
Π
n
×n
h
v
h
.
Π = IE
=
(xr,s ) (xvr,s )⊤
xvr,s
0nv ×nh
Πv
Let us now define the covariance of the state update as
("
#
#
)
"
h
i
h


x
Π
Π
hv
h
⊤
⊤
r+1,s
Π′ = IE
,
=
xhr+1,s
xvr,s+1
xvr,s+1
Π⊤
Πv
hv

(6)

(7)

⊤
⊤
nh ×nh
where Πhv = A1 Πh A⊤
,
3 + A2 Πv A4 + Qhv , Πvh = Πhv , and the dimensions are Πh ∈ R
nh ×nv
nv ×nh
nv ×nv
Πhv ∈ R
, Πvh ∈ R
, and Πv ∈ R
. The state covariance update equation
becomes

Π′ = AΠA⊤ + Q,

(8)

where Π = Π⊤ and Π′ = (Π′ )⊤ . Note that (8) is not a matrix Lyapunov state covariance
equation since Π′ 6= Π. However, by partitioning (8), one can decompose it into a pair of
coupled horizontal and vertical matrix Lyapunov type equations (Ramos & Mercère, 2016b).
Nevertheless, one can enforce the constraint Πhv = 0nh ×nv , which results in the joint matrix
Lyapunov equation





⊤ 

Πh
A1 A2
Πh
A1 A2
Qhh Qhv
0nh×nv
0nh ×nv
=
+
,
0nv ×nh
Πv
A3 A4
0nv ×nh
Πv
A3 A4
Qvh Qvv
or, more compactly,
Π = AΠA⊤ + Q.

(9)

Throughout the rest of this note we will use the symbol > 0 (≥ 0) to indicate that a matrix
is positive definite (positive semi-definite). Model (1a) – (1c) then satisfies the following
constraints, also known as the positive real conditions:
"
#
Q S
≥ 0, Q ≥ 0, R > 0, Π > 0.
(10)
S⊤ R
The 2-D output autocovariance sequence Λk,m ∈ Rny ×ny is given in terms of the Markov
parameters of the system as

⊤
⊤

if k = 0, m = 0

 C1 Πh C1 + C2 Πv C2 + R,


k−1

C1 A1 G1 ,
if k ≥ 1, m = 0
⊤
Λk,m = IE yr+k,s+m yr,s
=
(11)
m−1

G
,
if
k
=
0,
m
≥
1
C
A

2
2
4


 CAk−1,m G1,0 + CAk,m−1 G0,1 , if k ≥ 1, m ≥ 1,
where G1 and G2 are defined, respectively, as the horizontal and vertical partitions of the
matrix G ∈ Rnx ×ny , obtained from
)
#
("
xhr+1,s
⊤
= AΠC ⊤ + S,
(12a)
yr,s
G = IE
v
xr,s+1
2

with G1 ∈ Rnh ×ny and G2 ∈ Rnv ×ny given as
G1 = A1 Πh C1⊤ + A2 Πv C2⊤ + Sh
G2 = A3 Πh C1⊤ + A4 Πv C2⊤ + Sv ,
and
G1,0 =

"

G1
0nv ×ny

A0,0 = In ,

#

,

A1,0 =

G0,1 =
"

(12b)
(12c)
"

0nh×ny

A1

G2
A2

0nv ×nh 0nv ×nv

#

#

,

A0,1 =

"

0nh×nh 0nh ×nv
A3

A = A1,0 + A0,1 ,
Ak,m = A1,0 Ak−1,m + A0,1 Ak,m−1 ,
A−k,m = Ak,−m = 0nx ×nx , for k ≥ 1, m ≥ 1.

A4

#

for (k, m) > (0, 0)

The problem can now be stated as follows:
Definition 1. Given a data matrix Y ∈ Rny (N +1)×(M +1) corresponding to the output sequence
yr,s ∈ Rny , for r = 0, 1, . . . , N and s = 0, 1, . . . , M, find: (i) the system orders nh and nv
such that nx = nh + nv , (ii) parameter matrices {A, C, G} up to a similarity transformation,
v
N
(iii) covariance matrices {Π, Q, R, S}, and (iv) the initial conditions {xh0,s }M
s=0 and {xr,0 }r=0 ,
subject to the constraints (10), so that the 2nd-order statistics of the output of the system
match those of the given output data.
In order to simplify the analysis, we also formulate the problem in the innovations form
as
x
bhr+1,s = A1 x
bhr,s + A2 x
bvr,s + K1 er,s
x
bvr,s+1

=

yr,s =

A3 x
bhr,s
C1 x
bhr,s

(13a)

+ A4 x
bvr,s + K2 er,s
+ C2 x
bvr,s + er,s ,

(13b)
(13c)

where x
bhr,s ∈ Rnh and x
bvr,s ∈ Rnv are, respectively, the horizontal and vertical state estin
 o
h
h ⊤
∈ Rnh ×nh and Pv =
mates, with state estimate covariance matrices Ph = IE x
br,s x
br,s
n
o
n
⊤
⊤ o
IE x
bvr,s x
bvr,s
∈ Rnv ×nv . Furthermore, we assume that Phv = IE x
bhr,s x
bvr,s
= 0nh ×nv .
These state estimate covariance matrices satisfy the joint Riccati equation
P = AP A⊤ + (G − AP C ⊤ )(Λ0,0 − CP C ⊤ )−1 (G − AP C ⊤ )⊤ ,

where
P =

"

Ph

0nh ×nv

0nv ×nh

Pv

#

(14)

(15)

is 
a positive definite matrix. We further define the innovations covariance matrix Re =
IE er,s e⊤
r,s as
Re = Λ0,0 − C1 Ph C1⊤ − C2 Pv C2⊤

(16)

3

and state estimate errors and state estimate error covariance matrices, respectively, as
x̃hr,s = xhr,s − x
bhr,s ∈ Rnh
x̃vr,s = xvr,s − x
bvr,s ∈ Rnv
n
⊤ o
Σh = IE x̃hr,s x̃hr,s
= Πh − Ph ∈ Rnh ×nh
n
⊤ o
= Πv − Pv ∈ Rnv ×nv .
Σv = IE x̃vr,s x̃vr,s

Then Σh and Σv satisfy the joint Riccati equation

Σ = AΣA⊤ + Q + (AΣC ⊤ + S)(CΣC ⊤ + R)−1 (AΣC ⊤ + S)⊤ ,
where
Σ =

"

Σh

0nh ×nv

0nv ×nh

Σv

#

=

"

Πh

0nh ×nv

0nv ×nh

Πv

#

−

"

(17)

Ph

0nh ×nv

0nv ×nh

Pv

#

.

Finally, the Kalman gain matrix is given by either of the following two expressions
K = (G − AP C ⊤ )(Λ0,0 − CP C ⊤ )−1 ∈ Rnx ×ny
K = (AΣC ⊤ + S)(CΣC ⊤ + R)−1 ∈ Rnx ×ny ,
where
K =

"

K1
K2

#

(18a)
(18b)

,

with dimensions K1 ∈ Rnh ×ny and K2 ∈ Rnv ×ny .

2

Horizontal Data Processing

Let the horizontal and vertical past and future state matrices for k = 0, 1, . . . , M and N =
2i + j − 2 be defined as
 h

b h (k) ,
x
b0,k x
bh1,k x
bh2,k · · · x
bhj−1,k ∈ IRnh ×j
X
(19)
p


b h (k) ,
bhi+1,k x
bhi+2,k · · · x
bhi+j−1,k ∈ IRnh ×j
x
bhi,k x
(20)
X
f
 v

x
b0,k x
bv1,k x
bv2,k · · · x
bvj−1,k
v
 x
x
bvj,k 
bv3,k · · ·
bv2,k x
 bv1,k x

v
v
v

x
b
x
b
x
b
·
·
·
x
b
b vh (k) , 
2,k
3,k
4,k
j+1,k
(21)
X

 ∈ IRnv i×j
p
 .

.
.
.
.
..
..
..
 ..

..





vh
b
Xf (k) , 



x
bvi−1,k x
bvi,k x
bvi+1,k · · · x
bvi+j−2,k
x
bvi,k
x
bvi+1,k
x
bvi+2,k
..
.

x
bv2i−1,k

x
bvi+1,k
x
bvi+2,k
x
bvi+3,k
..
.
x
bv2i,k

x
bvi+2,k
x
bvi+3,k
x
bvi+4,k
..
.

x
bv2i+1,k

··· x
bvi+j−1,k
··· x
bvi+j,k
··· x
bvi+j+1,k
..
.
..
.
v
··· x
b2i+j−2,k
4






 ∈ IRnv i×j ,



(22)

where throughout the sequel, subscripts p and f denote past and future, respectively, superscripts h and v denote horizontal and vertical, respectively, vh denotes vertical from
horizontal data processing, and i and j are fixed integer constants such that j ≫ i and
ny i ≫ max{nh , nv }.
Likewise, we define the horizontal past and future innovations and output data matrices for
k = 0, 1, . . . , M and N = 2i + j − 2 as follows:


e0,k e1,k e2,k · · · ej−1,k
 e1,k e2,k e3,k · · ·
ej,k 




(23)
Eph (k) ,  e2,k e3,k e4,k · · · ej+1,k  ∈ IRny i×j

 ..
..
..
..
.
.

 .
.
.
.
.





Efh (k) , 





ei−1,k ei,k ei+1,k · · · ei+j−2,k
ei,k
ei+1,k
ei+2,k
..
.

ei+1,k
ei+2,k
ei+3,k
..
.

e2i−1,k

e2i,k

y0,k
 y1,k


h
Yp (k) ,  y2,k
 ..
 .
yi−1,k

yi,k
 yi+1,k

Yfh (k) , 
..

.

y2i−1,k

y1,k
y2,k
y3,k
..
.


· · · ei+j−1,k
· · · ei+j,k 

· · · ei+j+1,k 
 ∈ IRny i×j

.
.
..

..
e2i+1,k · · · e2i+j−2,k

y2,k · · · yj−1,k
y3,k · · ·
yj,k 

y4,k · · · yj+1,k 
 ∈ IRny i×j

.
..
.
..

..
.
ei+2,k
ei+3,k
ei+4,k
..
.

yi,k yi+1,k · · · yi+j−2,k
yi+1,k
yi+2,k
..
.
y2i,k

yi+2,k
yi+3,k
..
.

y2i+1,k

· · · yi+j−1,k
· · · yi+j,k
..
.
..
.
· · · y2i+j−2,k





 ∈ IRny i×j .


(24)

(25)

(26)

One can easily show that the following equations are satisfied for k = 0, 1, . . . , M
h h
bph (k) + Γvh
b vh
Yph (k) = Γhi X
i Xp (k) + Ki Ep (k)
h h
bfh (k) + Γvh
b vh
Yfh (k) = Γhi X
i Xf (k) + Ki Ef (k)

h h
vh b vh
b h (k) = Ai X
bh
X
f
1 p (k) + Φi Xp (k) + Li Ep (k),

h
where {Γhi , Φvh
i , Li } and other related matrices are defined as follows:


C1
 C1 A1 


h
Γi , 
 ∈ IRny i×nh
..


.

(28)
(29)

(30)

C1 Ai−1
1


i−2
Ai−1
∈ IRnh ×ny i
1 K1 A1 K1 · · · K1

 i−1
Ai−2
· · · Inh ∈ IRnh ×nh i
,
A1
1

Lhi ,

Φhi



(27)

5

(31)
(32)

Φvh
=
i
Θhi



i−2
Ai−1
1 A2 A1 A2 · · · A2


In h
 A1 


,  .  ∈ IRnh i×nh ,
.
 . 



∈ IRnh ×nv i = Φhi · (Ii ⊗ A2 )

(33)

(34)

Ai−1
1

and Ik denotes a (k×k) identity matrix. Finally, we define the lower triangular block Toeplitz
h
matrices {GhA1 , Γvh
i , Ki } as

GhA1





, 






Γvh
,

i


Kih



0ny ×nv
In h
0ny ×nv
..
..
..
.
.
.
i−2
i−3
A1
A1
· · · 0ny ×nv
C2
C1 A2
..
.

C2
..
.



 ∈ IRnh i×nh i


..

.
i−3
C1 Ai−2
A
C
A
A
·
·
· C2
2
1 1
2
1

= (Ii ⊗ C1 )GhA1 (Ii ⊗ A2 ) + (Ii ⊗ C2 )

In y
 C 1 K1
In y

, 
..
..
..

.
.
.

(35)





 ∈ IRny i×nv i


i−3
C1 Ai−2
1 K1 C1 A1 K1 · · · Iny

= (Ii ⊗ C1 )GhA1 (Ii ⊗ K1 ) + (Ii ⊗ Iny ).





 ∈ IRny i×ny i


(36)

(37)

For the purpose of horizontal data processing we will work with the equivalent horizontal
subsystem
x
bhr+1,s = A1 x
bhr,s + A2 x
bvr,s + K1 er,s
C1 x
bhr,s

yr,s =

+

C2 x
bvr,s

(38)

+ er,s ,

(39)

for r = 0, 1, . . . , N and s = 0, 1, . . . , M. However, at this point we need to make the following
notational simplification n̄h , nh (M + 1), n̄y , ny (M + 1), and ̄ , j(M + 1). Then, by
defining
Y hp ,
Y hf ,
ch
X
p

ch
X
f

cvh
X
p





Yph (0) Yph (1) · · · Yph (M)




∈ IRny ×̄

∈ IRny ×̄
i
h
b h (1) · · · X
b h (M) ∈ IRnh ×̄
b h (0) X
,
X
p
p
p
h
i
h
h
h
b
b
b
,
Xf (0) Xf (1) · · · Xf (M) ∈ IRnh ×̄
i
h
b vh (1) · · · X
b vh (M) ∈ IRnv i×̄
b vh (0) X
,
X
p
p
p
Yfh (0) Yfh (1) · · · Yfh (M)

6

(40)
(41)
(42)
(43)
(44)

i
b vh (1) · · · X
b vh (M) ∈ IRnv i×̄
b vh (0) X
X
f
f
f
 h

Ep (0) Eph (1) · · · Eph (M) ∈ IRny ×̄
,

 h
Ef (0) Efh (1) · · · Efh (M) ∈ IRny ×̄ ,
,

cvh ,
X
f
E hp

E hf

h

(45)
(46)
(47)

we get the horizontal subspace equations
h

vh

h

vh

ch
Ai1 X
p

cvh
Φvh
i Xp

c + Γvh X
c + K hE h
Y hp = Γhi X
p
i
p
i
p

(48)

c + Γvh X
c + K hEh
Y hf = Γhi X
f
i
f
i
f

2.1

ch
X
f

=

+

(49)

+ Lhi E hp .

(50)

Propagating the Vertical Hankel State Matrices

We will now propagate the state equation (1b) backward until we reach the initial vertical
states. By assuming zero initial vertical states, then the remaining vertical states are a
function of the innovations and horizontal states only. Since (21) and (22) are Hankel
matrices, we need to convert (1b) into a pair of past and future Hankel type matrix equations.
This is rather straightforward since (21) and (22) have partial horizontal dynamics (i.e., only
b h (k) and X
b h (k)). Thus, by substituting x
through X
bhr,s , x
bvr,s , and er,s in (1b) by their matrix
p
f
b h (k), X
b vh (k), E h (k)} and {X
b h (k), X
b vh (k), E h (k)}, we obtain, respectively,
equivalents, {X
p
p
p
f
f
f
the past and future vertical state equations given by
b vh (k + 1) = Θvh X
b h (k) + Avh X
b vh (k) + K vh E h (k)
X
p
i
p
i
p
i
p
vh
vh b h
vh b vh
vh h
b
Xf (k + 1) = Θi Xf (k) + Ai Xf (k) + Ki Ef (k),

(51)
(52)

where





Θvh
,

i

Avh
i

Kivh





, 


A3
A3 A1
..
.
A3 Ai−1
1
A4
A3 A2
..
.





 = (Ii ⊗ A3 ) · Θhi ∈ IRnv i×nh

A4
..
.



..

.
· · · A4
= (Ii ⊗ A3 )GA1 (Ii ⊗ A2 ) + (Ii ⊗ A4 )

K2
 A3 K1
K2

, 
..
..
..

.
.
.
A3 Ai−2
1 A2

(53)

A3 Ai−3
1 A2



 ∈ IRnv i×nv i


(54)



 ∈ IRnv i×ny i


(55)

i−3
A3 Ai−2
1 K1 A3 A1 K1 · · · K2
= (Ii ⊗ A3 )GA1 (Ii ⊗ K1 ) + (Ii ⊗ K2 ).

7



Let us now solve (51) and (52) recursively for k = 0, 1, . . . , M as follows:
b vh (0) = X
b vh (0)
X
p
p
vh
b
b
Xf (0) = Xfvh (0)
b vh (1) = Θvh X
b h (0) + Avh X
b vh (0) + K vh E h (0)
X
p

i

bfvh (1)
X

=

p

bh
Θvh
i Xf (0)

i

+

p

i

b vh
Avh
i Xf (0)

+

p

Kivh Efh (0)

b h (0) + (Avh )2 X
b vh (0) + K vh E h (1) + Avh K vh E h (0)
b vh (2) = Θvh X
b h (1) + Avh Θvh X
X
i
p
i
i
p
p
i
i
p
i
p
p
i
vh
vh
h
vh
vh
h
vh
2
vh
vh
h
vh
vh
b (2) = Θ X
b (1) + A Θ X
b (0) + (A ) X
b (0) + K E (1) + A K E h (0)
X
f
i
f
i
i
f
i
f
i
f
i
i
f
..
.
M
−1
M
−1
X
X
vh M −k−1 vh b h
M −k−1 vh h
b vh (M) = (Avh )M X
b vh (0) +
X
(A
)
Θ
X
(Avh
(k)
+
Ki Ep (k)
p
i
p
i
i
p
i )
b vh (M) = (Avh )M X
b vh (0) +
X
f
i
f

k=0
M
−1
X
k=0

M −k−1 vh b h
(Avh
Θi Xf (k) +
i )

(56)

k=0
M
−1
X

M −k−1 vh h
(Avh
Ki Ef (k). (57)
i )

k=0

b h (k), X
b h (k)} for k = 0, 1, . . . , M to construct upper triangular
Now we use {Eph (k), Efh (k), X
p
f
block Toeplitz matrices such as

 h
Eph (M)
Ep (0) Eph (1) · · ·

Eph (0) · · · Eph (M − 1) 


(58)
E ⋆p , 
 ∈ IRn̄y i×̄
..
.
.


.
.
E ⋆f

c⋆
X
p
c⋆
X
f



Eph (0)

Efh (M)
Efh (0) Efh (1) · · ·

Efh (0) · · · Efh (M − 1)

, 
..
..

.
.


Efh (0)

b h (0) X
b h (1) · · ·
b h (M)
X
X
p
p
p

h
h
b
b
Xp (0) · · · Xp (M − 1)

, 
..
..

.
.

h
b
Xp (0)

b h (0) X
b h (1) · · ·
b h (M)
X
X
f
f
f

h
h
b
b
Xf (0) · · · Xf (M − 1)

, 
..
..

.
.

h
b
X (0)
f





 ∈ IRn̄y i×̄





 ∈ IRn̄h ×̄






 ∈ IRn̄h ×̄ .



(59)

(60)

(61)

Notice that (58) – (59) contain block Hankel entries, thus are block Toeplitz with Hankel
blocks (BTHB). Finally, we define the controllability-like matrices
 vh

vh
M −1 vh
Θi Avh
· · · (Avh
Θi
Avh
∈ IRnv i×nh M
(62)
i Θi
i )
M ,
 vh

vh
nv i×ny iM
vh vh
vh M −1 vh
Ki
Ai Ki
· · · (Ai )
Ki
KM ,
∈ IR
(63)
8

∆Xbpvh (0) ,
∆Xb vh (0) ,
f

h

h

b vh (0) · · · (Avh )M X
b vh (0)
b vh (0) Avh X
X
i
p
i
p
p

b vh (0) Avh X
b vh (0) · · · (Avh )M X
b vh (0)
X
i
i
f
f
f

i
i

∈ IRnv i×̄

(64)

∈ IRnv i×̄ .

(65)

It can now be easily shown that the vertical states satisfy a pair of Hankel matrix equations
such as

 ⋆ 

cvh = ∆ b vh + 0nv i×n Avh X
c + 0nv i×ny i K vh E ⋆
X
(66)
p
p
p
M
M
h
Xp (0)
 ⋆ 


c + 0nv i×ny i K vh E ⋆ .
cvh = ∆ b vh + 0nv i×n Avh X
(67)
X
f
f
f
M
M
h
X (0)
f

b vh (0) = X
b vh (0) = 0nv i×j , then ∆ b vh = 0nv i×̄ and ∆ b vh = 0nv i×̄ .
If we now assume that X
p
f
Xp (0)
X (0)
f

cvh and X
cvh as
We then obtain the final expressions for X
p
f
cvh =
X
p

2.2

cvh
X
f

=





0nv i×nh Avh
M

0nv i×nh Avh
M




c⋆ +
X
p
c⋆
X
f

+





0nv i×ny i K vh
M

0nv i×ny i K vh
M





E ⋆p

(68)

E ⋆f .

(69)

Computing the Orthogonal Projection Y hf /Y hp

Since the vertical state matrices are now functions of horizontal states and innovations, we
can substitute these in the horizontal state equation (50), to get
c ⋆ = Ah X
c⋆ + K h E ⋆ ,
X
f
M
p
M p

where

AhM

K hM



vh
vh vh
Ai1 Φvh
Φvh
i Θi
i Ai Θi
i
vh

A1
Φi Θvh
i

i

A
1
, 





vh
vh vh
Lhi Φvh
Φvh
i Ki
i Ai Ki
h
vh

Li
Φi Kivh


Lhi
, 



(70)

vh M −1 vh
· · · Φvh
Θi
i (Ai )
vh
vh M −2 vh
· · · Φi (Ai )
Θi
vh M −3 vh
)
Θi
(A
· · · Φvh
i
i
.
..
..
.
Ai1
vh M −1 vh
· · · Φvh
Ki
i (Ai )
vh
vh M −2 vh
· · · Φi (Ai )
Ki
vh M −3 vh
· · · Φvh
Ki
i (Ai )
..
..
.
.
Lhi






 ∈ IRn̄h ×n̄h







 ∈ IRn̄h ×n̄y i .



(71)

(72)

Let us now re-visit (68) – (69) and further substitute (70) in (69), i.e.,
cvh = Q1 X
c ⋆ + Q2 E ⋆
X
p
p
p

(73)

cvh = Q1 X
c ⋆ + Q2 E ⋆
X
f
f
f


h c⋆
= Q1 AM X p + K hM E ⋆p + Q2 E ⋆f
c ⋆ + P 2 E ⋆ + Q2 E ⋆ ,
= P 1X
p
p
f

9

(74)

where
Q1 =
Q2 =
P1 =
P2 =











∈ IRnv i×nh M

∈ IRnv i×ny iM
K vh
M

· AhM ∈ IRnv i×n̄h
Avh
M

· K hM ∈ IRnv i×n̄y i .
Avh
M

0nv i×nh Avh
M
0nv i×ny i
0nv i×nh
0nv i×nh

(75)
(76)
(77)
(78)

We now compute the orthogonal projection Y hf /Y hp as
ch /Y h + Γvh X
cvh /Y h + K h E h /Y h
Y hf /Y hp = Γhi X
f
p
i
f
p
i
f
p


 

vh
1  h ch
⊤
⊤
c Y h + K h E h Y h ⊤ Rh −1 Y h ,
Γi X f Y hp + Γvh
X
=
i
f
p
i
f
p
pp
p
̄
where Rhpp = 1̄ Y hp Y hp

Y

h
h
f /Y p

⊤

. Furthermore, we now substitute Y hp

⊤

(79)

to get





−1 h
1 h ch ch ⊤ h ⊤ cvh ⊤ vh ⊤
h ⊤
h ⊤
Γi X f X p
Γi + X p
Γi
+ Ep
Ki
Rhpp
Yp
=
̄




−1 h
1 vh cvh ch ⊤ h ⊤ cvh ⊤ vh ⊤
h ⊤
h ⊤
+ Γi X f
Xp
Γi + X p
Γi
+ Ep
Ki
Rhpp
Yp
̄




−1 h
1 h h ch ⊤ h ⊤ cvh ⊤ vh ⊤
h ⊤
h ⊤
Γi + X p
Γi
+ Ep
Ki
Rhpp
Y p.
+ Ki E f X p
̄

Let us now look at each term individually. We start with

h
1 hc
Γ Xf
̄ i



ch
X
p

⊤

Γhi

⊤

  h ⊤
⊤
1 h ch ch ⊤ h ⊤
1 h  i ch
cvh + Lh E h X
c
Γi X f X p
Γi A1 X p + Φvh
X
Γi
=
Γhi
i
p
i
p
p
̄
̄
  h ⊤
⊤
1 h i ch ch ⊤ h ⊤ 1 h vh  c⋆
c
=
Γi A1 X p X p
Γhi
Γi + Γi Φi Q1 X p + Q2 E ⋆p X
p
̄
̄


⊤
⊤
1
ch
Γhi
+ Γhi Lhi E hp X
p
̄


Ph




 0nh ×nh  h ⊤
vh
h i
h ⊤
h vh
= Γi A1 Ph Γi + Γi Φi
0nv i×nh AM 
 Γi + 0ny i×ny i
..


.

h ⊤

= Γhi Ai1 Ph Γi

0nh ×nh

.

10

We continue with

h
1 hc
Γ
X
i
f
̄



cvh
X
p

⊤

Γvh
i

⊤

,

  vh ⊤

1 h ch cvh ⊤ vh ⊤
1 h  i ch
vh ⊤
cvh + Lh E h X
c
Γi X f X p
Γi A1 X p + Φvh
X
Γi
=
Γ
i
p
i
p
p
i
̄
̄
 

⊤


1 h i ch
c⋆ Q⊤ + E ⋆ ⊤ Q⊤ Γvh ⊤
=
Γi A1 X p
X
p
1
p
2
i
̄


⊤

1
vh ⊤
cvh X
cvh
X
Γ
+ Γhi Φvh
i
p
p
i
̄
 


⊤
⋆ ⊤
1 h h h
⊤
⋆ ⊤
⊤
c
+ Γi L i E p
X p Q1 + E p Q2
Γvh
.
i
̄
Thus we get,


i 0n ×n i
h
⊤
1 h ch cvh ⊤ vh ⊤
v
h
 Γvh
+ 0ny i×ny i
Γi X f X p
Γi
= Γhi Ai1 Ph 0nh ×nh · · · 0nh ×nh 

i
⊤
̄
Avh
M

vh ⊤
+Γhi Φvh
+ 0ny i×ny i
i (Ii ⊗ Pv ) Γi


i 0n i×n i
h
⊤
y
v
 Γvh
+Γhi Lhi Re 0ny ×ny · · · 0ny ×ny 

i
⊤
K vh
M

vh ⊤
= Γhi Φvh
.
i (Ii ⊗ Pv ) Γi
h

c Eh
Let us continue with the next term 1̄ Γhi X
f
p

⊤

Kih

⊤




⊤
⊤
1 h ch
1 h  i ch
h ⊤
vh cvh
h h
h ⊤
Γi X f E p
Γi A1 X p + Φi X p + Li E p E hp
Ki
=
Kih
̄
̄


⊤
⊤ 1
⋆
1 h i ch
⋆
c
Q
X
Γi A1 X p E hp
+
Q
E
Kih + Γhi Φvh
=
1
2 p
i
p
̄
̄
⊤
⊤
1
+ Γhi Lhi E hp E hp
Kih
̄

Re

  0ny ×ny

= 0ny i×ny i + 0ny i×ny i + Γhi Φvh
0nv i×ny i K vh

..
i
M

.

h ⊤

+Γhi Lhi (Ii ⊗ Re ) Ki
⊤
= Γhi Lhi (Ii ⊗ Re ) Kih .
11

0ny ×ny

E hp



⊤

Kih


⊤

 Kih


⊤


cvh Y h ⊤ , i.e.,
We then continue with the term 1̄ Γvh
i Xf
p





1 vh cvh
1 vh cvh ch ⊤ h ⊤ cvh ⊤ vh ⊤
h ⊤
h ⊤
h ⊤
Γ Xf Y p
Γ Xf
Xp
=
Γi + X p
Γi
+ Ep
Ki
̄ i
̄ i
⊤
⊤
1 vh cvh ch ⊤ h ⊤ 1 vh cvh cvh ⊤ vh ⊤ 1 vh cvh
Γi X f X p
Kih
Γi + Γi X f X p
Γi
+ Γi X f E hp
=
̄
̄
̄
  h ⊤

⊤

⊤
1 vh  c⋆
1
⊤
vh cvh cvh
h
c
=
Γi P 1 X p + P 2 E ⋆p + Q2 E ⋆f X
Γ
X
X
Γ
+
Γvh
i
f
p
i
p
i
̄
̄




1
c ⋆ + P 2 E ⋆ + Q2 E ⋆ E h ⊤ K h ⊤
+ Γvh
P 1X
i
p
p
f
p
i
̄


Ph

0nh ×nh 


 h ⊤
h 
vh
vh
= Γi
0nv i×nh AM · AM 
 Γi + 0ny i×ny i + 0ny i×ny i
..


.
0nh ×nh


⊤
⊤
1
cvh X
cvh
Γvh
+ 0ny i×ny i
+ Γvh
i Xf
p
i
̄


Re
 0ny ×ny 

⊤


h 
vh
·
K
+Γvh
0
A
 Kih + 0ny i×ny i
..
n
i×n
v
i
M 
M
h


.

0ny ×ny


1 vh cvh cvh ⊤ vh ⊤
=
Γ Xf Xp
Γi
.
̄ i
⊤
is zero since the future innovations are uncorrelated
Finally, the last term 1̄ Kih E hf Y hp
with the past data.
Now collecting all terms, we obtain

⊤

 
−1 h
vh ⊤
h
h ⊤
Y hf /Y hp = Γhi Ai1 Ph Γhi + Φvh
(I
⊗
P
)
Γ
+
L
(I
⊗
R
)
K
Rhpp
Yp
i
v
i
e
i
i
i
i
 vh ⊤
⊤
−1 h
1
cvh X
c
+ Γvh
X
Γvh
Rhpp
Y p.
i
f
p
i
̄

One can show that

⊤

⊤
vh ⊤
∆hi = Ai1 Ph Γhi + Φvh
+ Lhi (Ii ⊗ Re ) Kih
i (Ii ⊗ Pv ) Γi

 i−1
.
=
A1 G1 Ai−1
1 G1 · · · G1

Therefore, we have

Y hf /Y hp = Γhi · ∆hi Rhpp

−1

 vh ⊤
⊤
−1 h
1
cvh X
c
Y hp + Γvh
X
Γvh
Rhpp
Y p.
i
f
p
i
̄

We define the bias term as
−1 h
1 vh cvh cvh ⊤ vh ⊤
Γi X f X p
Γi
Rhpp
Y p.
bias =
̄
12




vh ⊤
c
Xp
We now need to find a closed form expression for
. For this, we will use
equations (29) and (51) – (52), along the following state estimate covariance equations
vh
1c
X
f
̄

Ph
Phv
Pvh
Pv

=
=
=
=

⊤
⊤
A1 Ph A⊤
1 + A2 Pv A2 + K1 Re K1
⊤
⊤
A1 Ph A⊤
3 + A2 Pv A4 + K1 Re K2
⊤
⊤
A3 Ph A⊤
1 + A2 Pv A2 + K2 Re K1
⊤
⊤
A3 Ph A⊤
3 + A4 Pv A4 + K2 Re K2 .

One can easily prove the following results.
h


 i
vh
i
vh ⊤
vh
vh ⊤
h
vh ⊤
P0 = Θi · A1 Ph Θi
+ Φ (Ii ⊗ Pv ) Ai
+ Li (Ii ⊗ Re ) Ki
h
= Θvh
i Φi (Ii ⊗ Phv )


⊤
vh ⊤
vh ⊤
Q = Θvh
+ Avh
+ Kih (Ii ⊗ Re ) Kivh
i Ph Θ i
i (Ii ⊗ Pv ) Ai


⊤
I
⊗
A
= (Ii ⊗ Pv ) + (Ii ⊗ A3 ) GA1 (Ii ⊗ Phv ) + (Ii ⊗ Pvh ) G⊤
i
3 .
A1
{z
}
|
Q0

 vh ⊤
cvh X
c
Now, since 1̄ X
can be represented as
f
p


⊤
1 X b vh
1 cvh cvh ⊤
vh
b
X (k) Xp (k)
X
Xp
=
̄ f
̄ k=0 f
M

Let us consider each product term for k = 0, 1, . . . , M. Starting with k = 0, we have that
1 b vh  b vh ⊤
X (0) Xp (0)
= 0nv i×nv i .
̄ f
Continuing with k = 1, we have

b vh (1) = Θvh X
b h (0) + Avh X
b vh (0) + K vh E h (0)
X
p
i
p
i
p
i
p
vh
vh b h
vh b vh
vh h
b
Xf (1) = Θi Xf (0) + Ai Xf (0) + Ki Ef (0)
b vh (0) + K vh E h (0).
b vh (0) + Θvh Lh E h (0) + Avh X
b h (0) + Θvh Φvh X
= Θvh Ai X
i

1

p

i

i

i

p

i

p


⊤
b vh (1) X
b vh (1) , we get
Now computing the covariance 1̄ X
p
f

i

f

i

f

 
⊤
⊤
1  vh b h
1 b vh  b vh ⊤
vh b vh
vh h
bph (0)
Xf (1) Xp (1)
Θi Xf (0) + Ai Xf (0) + Ki Ef (0)
X
Θvh
=
i
̄
̄


⊤
⊤
⊤
⊤
vh
vh
h
vh
b (0)
+ X
Ai
+ Ep (0)
Ki
p
h
 i


vh ⊤
h
vh ⊤
vh
vh ⊤
i
vh
+ Li (Ii ⊗ Re ) Ki
+ Φi (Ii ⊗ Pv ) Ai
= Θi A1 Ph Θi
h
= Θvh
i Φi (Ii ⊗ Phv )
= P0 .

13

Likewise, for k = 2, we have
b vh (2) = Θvh X
b h (1) + Avh Θvh X
b h (0) + (Avh )2 X
b vh (0) + K vh E h (1) + Avh K vh E h (0)
X
p
i
p
i
i
p
i
p
i
p
i
i
p
vh
vh
h
vh
vh
h
vh
2
vh
vh
h
vh
vh
b (2) = Θ X
b (1) + A Θ X
b (0) + (A ) X
b (0) + K E (1) + A K E h (0)
X
f
i
f
i
i
f
i
f
i
f
i
i
f
vh i b h
vh vh vh b h
vh vh vh b vh
vh vh vh h
= Θi A1 Xp (1) + Θi Φi Θi Xp (0) + Θi Φi Ai Xp (0) + Θi Φi Ki Ep (0)
h h
vh vh i b h
vh vh vh b vh
vh vh h h
+Θvh
i Li Ep (0) + Ai Θi A1 Xp (0) + Ai Θi Φi Xp (0) + Ai Θi Li Ep (0)
2 vh
b (0) + Avh K vh E h (0) + K vh E h (1)
+ Avh
X
i
f
i
i
f
i
f


⊤
b vh (2) X
b vh (2) , we get
Now computing the covariance 1̄ X
p
f

h


 i
1 b vh  b vh ⊤
i
vh ⊤
vh
vh ⊤
h
vh ⊤
Xf (2) Xp (2)
= Θvh
A
P
Θ
+
Φ
Q
A
+
L
(I
⊗
R
)
K
h
i
e
i
1
i
i
i
i
̄

vh h
vh ⊤
+Avh
i Θi Φi (Ii ⊗ Phv ) Ai
h


 i
i
vh ⊤
vh
vh ⊤
h
vh ⊤
A
P
+
L
(I
⊗
R
)
K
= Θvh
Θ
+
Φ
(I
⊗
P
)
A
i
v
i
e
1 h
i
i
i
i
i


vh
vh ⊤
vh h
vh ⊤
+Θvh
+ Avh
i Φi Q0 Ai
i Θi Φi (Ii ⊗ Phv ) Ai


h
vh vh
vh ⊤
vh h
vh ⊤
= Θvh
+ Avh
i Φi (Ii ⊗ Phv ) + Θi Φi Q0 Ai
i Θi Φi (Ii ⊗ Phv ) Ai


vh ⊤
vh vh
vh ⊤
= P0 + Avh
P
A
+
Θ
Φ
Q
A
.
0
0
i
i
i
i
i

Continuing further, for k ≥ 3, we obtain the general expression for
as

vh
1 vh c
Γ
X
i
f
̄

 vh ⊤
⊤
c
X
Γvh
p
i

M
−ℓ X
M

X

 ⊤ vh ⊤
1 vh cvh cvh ⊤ vh ⊤
vh k
vh k
Γi
(M − ℓ − k + 1)Γvh
A
P
A
Γi X f X p
Γi
=
0
i
i
i
̄
k=0 ℓ=1

+

MX
−ℓ−1 M
−1
X

(M − ℓ − k)Γvh
Avh
i
i

k=0

×Q0



ℓ=1

Avh
i

ℓ+k ⊤

Γvh
i

⊤

.

ℓ−1

vh
Θvh
Avh
i Φi
i

k

(80)

Analyzing (80) one can see that each term in the first sum is a function of P0 and each term
in the second sum is a function of Q0 , both of which are functions of Phv and/or Pvh , which
by (15) are zero matrices. Thus, we conclude that the bias term is zero. Thus,

and

−1 h
1 vh cvh cvh ⊤ vh ⊤
bias =
Γi X f X p
Γi
Rhpp
Y p = 0ny i×jh
̄
Y hf /Y hp = Γhi · ∆hi Rhpp
h

c .
= Γhi · X
f

−1

Y hp
14

2.3

Improving The State Estimates

Since the orthogonal projection is not exact, there is a small bias introduced that may affect
the identification of the system parameters. Despite the fact that the bias is rather small,
one can iterate the procedure in order to improve the state estimates and eliminate the bias.
We now propose an oblique projection approach to improve the state estimates. Along the
way we also propose a procedure for computing the initial states. We start by assuming that
the vertical states1 are available from an orthogonal projection in the vertical direction, i.e.,
cv (see Ramos and Mercère (2017a) for details), where X
cv ∈ Rnv ×̄ . Then
Y vf /Y vp ∼
= Γvi · X
f
f
cvh ∈ Rnv i×̄ .
we assemble the vertical from horizontal data processing Hankel state matrix X
f

The second stage of the algorithm starts by defining W hp as
 vh 
c
X
p
 ∈ IR(nv +ny )i×̄ .
W hp = 
Y hp

We then compute the RQ decomposition of the past/future data as follows:
 vh 



⊤
c
X
Q
R
 f 
 11
 1 

 ⊤ 
h  = 
 Wp 
 R21 R22
  Q2  ,





h
⊤
Yf
R31 R32 R33
Q3

(81)

(82)

where R11 ∈ Rnv i×nv i , R21 ∈ R(nv +ny )i×nv i , R22 ∈ R(nv +ny )i×(nv +ny )i , R31 ∈ Rny i×nv i , R32 ∈
Rny i×(nv +ny )i , R33 ∈ Rny i×ny i , Q1 ∈ Rj×nv i , Q2 ∈ Rj×(nv +ny )i , and Q3 ∈ Rj×ny i .
From (86) one can find an expression for Y hf using the R and Q parameters, along with

cvh and Q⊤ = R−1 W h − R21 Q⊤ . That is,
using Q⊤ = R−1 X
1

11

f

2

22

1

p

⊤
⊤
Y hf = R31 Q⊤
1 + R32 Q2 + R33 Q3


h
−1
⊤
⊤
+
R
R
W
−
R
Q
= R31 Q⊤
32
21
1
22
p
1 + R33 Q3
 −1
−1
−1
cvh + R33 Q⊤
= R32 R22
W hp + R31 − R32 R22
R21 R11
·X
f
3
ch + Γvh · X
cvh + K h · E h .
= Γhi · X
f
i
f
i
f

It is now clearly evident that

ch = R32 R−1 W h
Γhi · X
f
22
p
vh

c
Γvh
i · Xf

=

 −1
−1
cvh
R31 − R32 R22
R21 R11
·X
f

Kih · E hf = R33 Q⊤
3.

h
Without computing the system parameters, our aim here is to compute Γhi , then Γvh
i and Ki
h
with the right lower triangular Toeplitz structure. Computing Γi is straight forward, thus
h
we assume it is already known. We now concentrate on computing Γvh
i and Ki .
1

Here we assume that one can compute the entire vertical state sequence x
bvr,s , for r = 0, 1, . . . , N and
s = 0, 1, . . . , M .

15

Using the Lower Triangular Toeplitz System Solver (LTTSS) procedure in Ramos and Mercère
(2016b), we compute Γvh
i by solving the linear system of equations
Iny i · Γvh
i · R11 =


−1
R31 − R32 R22
R21 ,

subject to Γvh
i being lower triangular Toeplitz. The solution is

−1
Γvh
= LTTSS{Iny i , R11 , R31 − R32 R22
R21 , ny , nv , i}.
i

We now define E f and E f1 as



E
(0)
E
(1)
·
·
·
E
(M)
E f = R33 Q⊤
=
∈ Rny i×̄
f
f
f
3


Ef1 (0) Ef1 (1) · · · Ef1 (M) ∈ Rny i×(j−i+1)(M +1) ,
E f1 =

where






Ef (k) = 


=



ek0,0
ek1,0
ek2,0
..
.

ek0,1
ek1,1
ek2,1
..
.

ek0,2
ek1,2
ek2,2
..
.

ek0,j−i+1
ek1,j−i+1
ek2,j−i+1
..
.

ek0,j−i
ek1,j−i
ek2,j−i
..
.

···
···
···
..
.

···
···
···
..
.

ek0,j−1
ek1,j−1
ek2,j−1
..
.

eki−1,0 eki−1,1 eki−1,2 · · · eki−1,j−i eki−1,j−i+1 · · · eki−1,j−1

Ef1 (k) ×









and × denotes a matrix that is not relevant to the discussion. Furthermore, since the main
diagonal blocks of Kih are all equal to Iny and all elements above the main diagonal blocks
are 0ny ×ny , we observe that the first ny rows of E f contains a sequence of innovations, from
which Kih can be computed. That is,


Iny 0ny ×ny · · · 0ny ×ny

where



e0 (k) =



Ef =

ek0,0 ek0,1 · · · ek0,j−1




∼
=

e0 (0) e0 (1) · · · e0 (M)



ei,k ei+1,k · · · ei+j−1,k




∈ Rny ×̄ ,

∈ Rny ×j .

Let us now form the array of Hankel matrices using e0 (k), for k = 0, 1, . . . , M, i.e.,
E f2 =
where



Ef2 (0) Ef2 (1) · · · Ef2 (M)





Ef2 (k) = 



ek0,0
ek0,1
ek0,2
..
.

ek0,1
ek0,2
ek0,3
..
.

ek0,2
ek0,3
ek0,4
..
.

ek0,i−1 ek0,i ek0,i+1



∈ Rny i×(j−i+1)(M +1) ,

· · · ek0,j−i
· · · ek0,j−i+1
· · · ek0,j−i+2
..
.
..
.
k
· · · e0,j−1
16






 ∈ Rny i×(j−i+1) .



Notice that if we knew E hf , then Ef2 (k) would be the first j − i + 1 columns of Efh (k), for
k = 0, 1, . . . , M. That is,


ej+1,k · · · ei+j−1,k
ei,k
ei+1,k ei+2,k · · ·
ej,k
 ei+1,k ei+2,k ei+3,k · · · ej+1,k ej+2,k · · · ei+j,k 




Efh (k) =  ei+2,k ei+3,k ei+4,k · · · ej+2,k ej+3,k · · · ei+j+1,k 


.
..
.
..
.
..
.
.
..
..
..


..
..
.
.
.
e2i−1,k e2i,k e2i+1,k · · · ei+j−1,k ei+j,k · · · e2i+j−2,k


Ef2 (k) × .
=

Let us now define the covariance matrices V 1 ∈ Rny i×ny i and V 2 ∈ Rny i×ny i as
1
E f (E f2 )⊤
(j − i + 1)(M + 1) 1
1
E f (E f2 )⊤ .
=
(j − i + 1)(M + 1) 2

V1 =
V2

We can now find a relationship between V 1 and V 2 as follows:
V 1 = Iny · Kih · V 2 ,
and upon applying the vec operator on both sides, we obtain

vec{V 1 } = V T2 ⊗ Iny i · vec{Kih }.

(83)

Now, since Kih is a lower triangular Toeplitz matrix, vec{Kih } will contain repeated elements.
To remove these redundancies, we apply the identity
vec{Kih } = FKih · kih ,
2 2

2

where FKih ∈ Rny i ×ny i is a permutation matrix with elements equal to 0 and 1 and kih
contains all the elements of Kih , i.e.,


vec{k0 }
 vec{k1 } 


h
ki = 
,
..


.
vec{ki−1 }

with

ks =



if s = 0
In y ,
C1 A1s−1 K1 , if s ≥ 1.

The (LTTSS) procedure performs this operation and assembles the full Toeplitz matrix.
Thus, we obtain Kih as
V = LTTSS{Iny i , V 2 , V 1 , ny , ny , i}

Kih = V Ii ⊗ K0−1 ,

17

where K0 ∈ Rny ×ny is the first (ny × ny ) block of V . Knowing Kih , we can now compute E hf
from R33 Q⊤
3 . For k = 0, 1, . . . , M, we need to solve for Ef (k) using the Hankel System Solver
(HSS) procedure outlined in Ramos and Mercère (2016a). That is, we solve the following
linear system of equations
Ef (k) = Kih · Efh (k) · Ij , for k = 0, 1, . . . , M,
subject to Ef (k) being a Hankel matrix. Upon applying the vec operator on both sides, we
get

vec{Ef (k)} = Ij ⊗ Kih vec{Efh (k)}.
(84)
However, since Efh (k) is a Hankel matrix, vec{Efh (k)}, will be inefficient for solving (84)
because of the repeated elements. In order to compute the minimum number of elements
from vec{Efh (k)}, we use the property
vec{Efh (k)} = FEfh(k) · ehf (k),
where FEfh(k) ∈ Rny ij×ny (i+j−1) is

vec{ei,k }
 vec{ei+1,k }

ehf (k) = 
..

.

a permutation matrix with elements equal to 0 and 1 and


vec{e2i+j−2,k }



.


The HSS procedure handles the removal of repeated elements and assembles the full Hankel
matrix EFh (k). Thus, we get
Efh (k) = HSS{Kih , Ij , Ef (k), ny , 1, i, j}, for k = 0, 1, . . . , M.
Once we have all M + 1 solutions, we can then assemble the full matrix E hf .
cvh using {Y h , X
ch , E h , Γh , Γvh , K h }. That is, let us define
We now need to find X
f
f
f
f
i
i
i
ch − K h E h = Γvh X
b vh = Y h − Γh X
cvh
Z
f
f
i
f
i
f
i
f
i
h
=
Zbfvh (0) Zbfvh (1) · · · Zbfvh (M) ∈ Rny i×̄ .

Once again, we have M + 1 systems of equations of the form
b vh
Zbfvh (k) = Γvh
i · Xf (k) · Ij , for k = 0, 1, . . . , M.

If we now apply the vec operator on both sides, we get

bvh
vec{Zbfvh (k)} = Ij ⊗ Γvh
· FXb vh (k) · x
f (k),
i
f

where FXb vh (k) ∈ Rnv ij×nv (i+j−1) is
f

vec{b
xvi,k }
 vec{b
xvi+1,k }

(k)
=
x
bvh

..
f

.

(85)

a permutation matrix with elements equal to 0 and 1 and


vec{b
xv2i+j−2,k }



.


18

By applying the (HSS) procedure, we obtain
b vh (k) = HSS{Γvh , Ij , Zbvh (k), nv , 1, i, j}, for k = 0, 1, . . . , M.
X
f
i
f

In horizontal data processing, as it relates to future data, we have two instances where
we need to use the HSS procedure. This operation could be computationally expensive since
it has to be done (M + 1) times. Nevertheless, the solution will give us the right structure
cvh . Now that we have X
b vh (k) for k = 0, 1, . . . , M, we can assemble the full
for E h and X
f

f

f

cvh .
X
f

matrix
Let us now define the following matrices:
h
 i
†
 h
h †
h †
∈ Rnh ×(nv +2ny )i
T h2 =
Γ
− Γhi Γvh
−
Γ
K
i
i
i
i
 vh 
c
X
 f 


H hf =  E hf  ∈ R(nv +2ny )i×̄ ,


h
Yf

†
ch
where Γhi is the pseudo-inverse of Γhi . Notice that the future horizontal state matrix X
f
is related to T h2 and H hf via
ch = T h H h .
X
f
2
f

The next step is
 vh 
c
X
f
 vh 
 X
c 
 p  =


 Y hp 
Y hf

ch . Let us re-visit the QR decomposition of the data, i.e.,
to recover X
p



R11
Q⊤
1
 1
 ⊤ 
1
 R21 R22
  Q21 



(86)
 R2 R2 R3
  Q⊤  ,
 21
  22 
22
22
1
2
R31 R32
R32
R33
Q⊤
3

1
2
1
2
3
where R21
∈ Rnv i×nv i , R21
∈ Rny i×nv i , R22
∈ Rnv i×nv i , R22
∈ Rny i×nv i , R22
∈ Rny i×ny i ,
1
ny i×nv i
2
ny i×ny i
j×nv i
j×ny i
⊤
R32 ∈ R
, R32 ∈ R
, Q21 ∈ R
, and Q22 ∈ R
. Then Q21 and Y hp can be
expressed as


1 −1 cvh
1 −1 1
Q⊤
=
R
X
−
R
R21 Q⊤
21
22
p
22
1
2
2
⊤
3
⊤
Y hp = R21
Q⊤
1 + R22 Q21 + R22 Q22




2
2
1 −1 cvh
1 −1 1
⊤
3
= R21
Q⊤
+
R
R
X
−
R
R
Q
+ R22
Q⊤
1
22
22
p
22
21 1
22




2
1 −1 cvh
3
2
2
1 −1 1
X p + R22
Q⊤
= R21
− R22
R22
R21 Q⊤
1 + R22 R22
22 .

vh
2
1 −1 c
We can now isolate R22
(R22
) X p from the rest. That is, we define




vh
h
3
⊤
2
1 −1 cvh
2
2
1 −1 1
b
R21 Q⊤
Z p = Y p − R21 − R22 R22
·X
1 − R22 Q22 = R22 R22
|
{z
} p

=

h

Zbpvh (0) Zbpvh (1) · · · Zbpvh (M)

i

19

= Γvh
i ·

h

Γvh
i

b vh (0) X
b vh (1) · · · X
b vh (M)
X
p
p
p

i

,

b vh ∈ Rny i×̄ and Zbvh (k) ∈ Rny i×j . It is clear that we have M + 1 equations of the
where Z
p
p
form
b vh
Zbpvh (k) = Γvh
i · Xp (k) · Ij , for k = 0, 1, . . . , M.

If we apply the vec operator on both sides of (87), we get

vec{Zbpvh (k)} = Ij ⊗ Γvh
FXbpvh (k) · x
bvh
i
p ,

(87)

(88)

where FXbpvh (k) ∈ Rnv ij×nv (i+j−1) is a permutation matrix with elements equal to 0 and 1 and


x
bv0,k


v
 x

b
1,k


x
bvh
(k)
=

.
..
p


.


v
x
bi+j−2,k
Now applying the HSS procedure, which solves (88) and reconstructs the Hankel matrix, we
can find the individual solutions from
b vh (k) = HSS{Γvh , Ij , Zbvh (k), nv , 1, i, j}, for k = 0, 1, . . . , M.
X
p
i
p

cvh matrix.
b vh (k), for k = 0, 1, . . . , M, we can now assemble the full X
Knowing all X
p
p
We now let



⊥  h
cvh = Γh ⊥ Γh X
c h + Γh ⊥ K h E h
·
X
Y p − Γvh
Y p = Γhi
p
i
i
i
i
i
p
p
{z
}
|


0(ny i−nh )×̄


⊥

Yp (0) Yp (1) · · · Yp (M) = Γhi Kih Eph (0) Eph (1) · · · Eph (M) ,
⊥
where Yp (k) ∈ R(ny i−nh )×j and Γhi
∈ R(ny i−nh )×ny i is the orthogonal complement of Γhi .
We now have M + 1 equations of the form
⊥
Yp (k) = Γhi Kih · Eph (k) · Ij , for k = 0, 1, . . . , M.
(89)
=



Let us now apply the vec operator on both sides of (89), to get


⊥
vec{Yp (k)} = Ij ⊗ Γhi Kih · FEph(k) · ehp (k), for k = 0, 1, . . . , M,

where FEph(k) ∈ Rny ij×ny (i+j−1) is a permutation matrix with elements equal to 0 and 1 and


vec{e0,k }
 vec{e1,k } 


ehp (k) = 
.
..


.
vec{ei+j−2,k }

Once again, the solution can be found by applying the HSS procedure,
⊥
Eph (k) = HSS{ Γhi Kih , Ij , Yp (k), ny , 1, i, j}, for k = 0, 1, . . . , M,
20

from which one can assemble the full E hp matrix.
cvh and E h , one can now compute X
ch from
Knowing X
p
p
p
ch =
X
p
=

Γhi

h


†  h
cvh − K h E h
Y p − Γvh
X
i
p
i
p

−

†
Γhi

Γvh
i

−

†
Γhi

Kih

†
Γhi

= T h2 H hp ,



cvh
X
p

i

 E hp

Y hp







where

H hp



cvh
X
p



=  E hp

Y hp





 ∈ R(nv +2ny )i×̄ .


Let us now define J ∈ Rnh ×(nh +(nv +ny )i) and H ∈ R(nh +(nv +ny )i)×̄ as


h
c
Xp


i
h

h
cvh 
and
H
=
J =
.

L
X
Ai1 Φvh
i
i
 p 
E hp
ch = J H, solving for J would require solving a large system of equations. In
Then since X
f
order to avoid this, we compute the covariances Z 1 and Z 2 as
1
HH ⊤ ∈ R(nh +(nv +ny )i)×(nh +(nv +ny )i)
̄
1 ch ⊤
=
X H ∈ Rnh ×(nh +(nv +ny )i) .
̄ f

Z1 =
Z2

Then, the solution for J becomes
J = Z 2 Z −1
1 ,

h
from which Ai1 , Φvh
i , and Li can be computed, i.e.,

Ai1 = J (:, 1 : nh )
Φvh
= J (:, nh + 1 : nh + nv i)
i
h
Li = J (:, nh + nv i + 1 : nh + (nv + ny )i).
21

ch into X
ch we get
Now substituting X
p
f
ch = Ai X
ch + Φvh X
cvh + Lh E h
X
f
1
p
i
p
i
p
=

=

Ai1

h

h

−

†
Γhi

Γvh
i

i
h
Φvh
i − A1 Γi

†

−

†
Γhi

†
Γhi

Kih

Γvh
Lhi − Ai1 Γhi
i

†



cvh
X
p

i

 E hp

Y hp





h h
cvh
 + Φvh
i X p + Li E p


Kih Ai1 Γhi

= T h1 H hp ,



cvh
X
p

† i 

 E hp

Y hp







where T h1 ∈ Rnh ×(nv +2ny )i is defined as
h
† vh
† h
† i
h
vh
i
h
h
i
h
i
h
T1 =
.
Φi − A1 Γi Γi Li − A1 Γi Ki A1 Γi
h

c from
The final task involves computing X
f+
h
i
h
h
h
n ×̄
c
h
h
h
b
b
b
X
=
T
H
=
Xf + (0) Xf + (1) · · · Xf + (M) ∈ R h ,
f+
1
f
where

b h (k) =
X
f+

h

bh2i+1,k · · · x
bh2i+j−1,k
x
bh2i,k x

i

∈ Rnh ×j .

b h (k) and X
b h (k) there is an overlap of {b
bhj−1,k }. That is,
bh2i+1,k , . . . , x
Now, between X
xh2i,k , x
p
f+


h
h
h
h
h
h
h
h
h
h
x
b
x
b
·
·
·
x
b
x
b
x
b
x
b
··· x
b2i−1,k x
b2i,k x
b2i+1,k · · · x
bj−1,k
j,k
j+1,k
2i+j−1,k
 | 0,k 1,k 2,k

{z
}
b h (k)
X
p

and


bh1,k x
bh2,k · · · x
bh2i−1,k x
x
bh0,k x
bh2i,k x
bh2i+1,k · · · x
bhj−1,k x
bhj,k x
bhj+1,k · · · x
bh2i+j−1,k

|
{z
} .
b h (k)
X
f+

Therefore, the entire horizontal state sequence can be recovered from
h
i
b h (:, 0 : 2i − 1)(k) X
b h (k) , for k = 0, 1, . . . , M.
X
p

f+

h

c by vectorizing the
Now form the (nh (N + 1) × (M + 1)) matrix of horizontal states X
h
h
b (k) and X
b (k), for k = 0, 1, . . . , M, then recalling that N = 2i+j −2,
individual terms in X
p
f+
i.e.,
i
h
n j×(M +1)
h
b
h
h
h
b
b
b
(90)
Xp ,
vec{Xp (0)} vec{Xp (1)} · · · vec{Xp (M)} ∈ R h
h
i
b h (0)} vec{X
b h (1)} · · · vec{X
b h (M)} ∈ Rnh j×(M +1) ,
Xbfh+ ,
(91)
vec{X
f+
f+
f+
22

ch from
and finally compute X


Xbph (0 : 2i − 1, :)
h
c , 

X
h
b
Xf + (0 : j − 2, :)

x
bh0,0
x
bh0,1 · · · x
bh0,M


h
 x
x
bh1,1 · · · x
bh1,M
 b1,0

..
..
..
..

.
.
.
.


h
 x
bh2i−1,1 · · · x
bh2i−1,M
 b2i−1,0 x
= 
 x
h
x
bh2i,1 · · · x
bh2i,M
 b2i,0

 x
h
bh2i+1,1 · · · x
bh2i+1,M
 b2i+1,0 x

..
..
..
..

.

.
.
.

x
bhN,0
x
bhN,1 · · · x
bhN,M













 ∈ Rnh (N +1)×(M +1) .










(92)

This completes the computation of the horizontal states. The same procedure must be
applied in the vertical direction to get the vertical state estimates.

References
Ramos, J. A., & Mercère, G. (2016a, October). Image modeling based on a 2-D stochastic subspace system
identification algorithm. Multidimensional Systems and Signal Processing, 28 , 1133–1165.
Ramos, J. A., & Mercère, G. (2016b). Subspace algorithms for identifying separable-in-denominator 2d
systems with deterministic-stochastic inputs. International Journal of Control , 89 (12), 2584–2610.
Retrieved from http://dx.doi.org/10.1080/00207179.2016.1172258
Ramos, J. A., & Mercère, G. (2017a). A stochastic subspace system identification algorithm for state
space systems in the general 2-d roesser model form. International Journal of Control (submitted for
publication).

23

