1

QUARKS: Identification of large-scale Kronecker
Vector-AutoRegressive models

arXiv:1609.07518v2 [] 20 Nov 2017

Baptiste Sinquin1 and Michel Verhaegen1

Abstract—In this paper we address the identification of 2D
spatial-temporal dynamical systems described by the VectorAuto-Regressive (VAR) form. The coefficient-matrices of the VAR
model are parametrized as sums of Kronecker products. When
the number of terms in the sum is small compared to the
size of the matrices, such a Kronecker representation efficiently
models large-scale VAR models. Estimating the coefficient matrices in least-squares sense gives rise to a bilinear estimation
problem which is tackled using an Alternating Least Squares
(ALS) algorithm. Regularization or parameter constraints on the
coefficient-matrices allows to induce temporal network properties
such as stability, as well as spatial properties such as sparsity or
Toeplitz structure. Convergence of the regularized ALS is proved
using fixed-point theory. A numerical example demonstrates the
advantages of the new modeling paradigm. It leads to comparable
variance of the prediction error with the unstructured leastsquares estimation of VAR models. However, the number of
parameters grows only linearly with respect to the number of
nodes in the 2D sensor network instead of quadratically in the
case of fully unstructured coefficient-matrices.
Index Terms—system identification, Vector Auto-Regressive
model, large-scale networks, Kronecker product, Alternating
Least Squares.

I. I NTRODUCTION
System identification of large-scale networks has received
an increased attention during the past years. Major scientific
and engineering projects such as [1] challenge the identification community to derive scalable identification algorithms.
In this work, the focus is on regular 2D sensor networks
defined on a square or circular rastered grid of size N × N
with N large. Such networks occur in multi-dimensional signal
processing problems, such as image/video processing [2] but
also in control for high resolution imaging using adaptive
optics [3]. Prior knowledge on the network structure is often of
prime importance to cope with the challenges inherent to the
large dimensions. Potential model structures for representing
the spatial-temporal dynamics include 2D state-space models,
e.g [4] [5], although the identification of the latter model
in a system identification or convex optimization context is
still an open research question. The difficulties that arise in
e.g the Roesser model is that spatial causality is not present,
hence making it a challenge to derive efficient identification
methods and global convergence under generic conditions. For
specific conditions, such as assuming a separability condition
1 Both authors are with the Delft Center for Systems and Control,
Technische Universiteit Delft, the Netherlands.
Corresponding author: baptiste.sinquin@gmail.com
The research leading to these results has received funding from the
European Research Council under the European Union’s Seventh Framework
Programme (FP7/2007-2013) / ERC grant agreement 339681.

of the transfer function called Causal, Recursive, Separable in
Denominator (CRSD), a subspace algorithm is derived in [6].
The temporal dimension as well as considering spatial varying
dynamics in a global (modal) context is not investigated
however. The global context is of key relevance to impose
constraints such as stability of the global network. Another
framework consists in assuming that each node is connected
to very few other nodes in the network. Identification of
these sparse networks has been investigated in [7] within
the Bayesian framework. Sequentially Semi-Separable (SSS)
networks are composed of a 1D string of interconnected
systems and have been analyzed in [8]. Each of the subsystem
is modeled with a mixed causal anti-causal linear time varying
model and shares unknown interconnections with the closest
neighbors. This framework reveals to be very elegant to
deal in a scalable manner with large strings of subsystems:
both linear algebra operations and control to achieve global
H2 performance were shown to be achievable within linear
computational complexity in the string’s size. Subspace system
identification for such systems has been achieved locally in
[9]. The extension of the 1D SSS methods to higher spatial
dimensions gives rise to multi-level SSS problems, for which
up till now no efficient solution for identification (and control)
exist. An alternative is proposed in [10] that introduces an
identification method using 2D Finite Impulse Response models. However, this identification and related methods search
for a local model estimation and have difficulties in assuming
and/or imposing global network properties such as the stability
of the overall system. Rather than having a zonal representation of networks as in the previous paragraph, the work in
[11] sheds the light onto the modal subsystems associated
with the network. A generalization is found in [12] which
introduces α-decomposable networks that allow for α different
kinds of subsystems in the whole network to interact. In
this paper, we present a modal network representation that
includes α-decomposable systems as a special case and allow
the subsystems to be heterogeneous.
In this paper, a novel modeling and identification paradigm
is introduced to model 2D spatial systems with temporal
dynamics. As a fundament of this new approach, we restrict
to temporal Vector Auto-Regressive (VAR) models with the
spatial structure imposed on the coefficient matrices. Let a
2D network connected on a grid of size N × N . The spatial
structure represents the coefficient matrices {Ai }i=1..p as a
finite sum of a Kronecker product between low dimensional
matrices:
r
X
2
2
Ai =
Ui ⊗ Vi ∈ RN ×N
i=1

2

where r  N is called the Kronecker rank and Ui , Vi ∈
RN ×N are the factor matrices. Such representation of large
dimension matrices was studied in [13] in which the equivalence between expressing a matrix as a sum containing
few Kronecker products and a low-rank approximation of a
reshuffled matrix was established. More than only enjoying
the storage of a reduced number of entries, such a structure
enables fast computations thanks to the very pleasant algebra
of the Kronecker product, see e.g [14].
Using Kronecker structures for forward modeling of multidimensional problems is well-known, especially in tensorbased scientific computing [15]. Besides, Kronecker structures have been applied efficiently for computing second
moments in multi-dimensional processes [17], for analyzing
EEG signals [16] and for image deblurring, [18]. The latter
example enables to relate the Kronecker rank-one modeling
with physical properties of the system. Denoting an object
O imaged with a static optical system, the resulting blurred
image B undergoes the linear blurring operation as follows:
vec(B) = Avec(O)

(1)

The coefficient-matrix A is related to the Point-Spread Function (PSF) (or 2D impulse response) of the optical system. The
equation (1) represents the 2D convolution operation between
the PSF and the object O. The structure in A is related to the
separability of the PSF, which implies the following Kronecker
structure for the coefficient-matrix A:
A = Ar ⊗ Ac

(2)

where Ar and Ac represent respectively the 1D convolution
with the rows and columns. A large-scale static input-output
map in (1) is represented by a Kronecker matrix as in (2). In a
more general context, separation-of-variable techniques have
been applied in [19] and the references therein to break down
the curse of dimensionality when modeling high-dimensional
partial differential equations.
Although tensor-based algorithms for handling large
datasets receive a growing interest, system identification of
multi-dimensional systems is however in its infancy. An
overview of data-driven algorithms that handle large datasets
using the tensor representation was provided in [22] among
which a multilinear tensor regression for relational longitudinal data, [23]. The approach proposed in [23] handles
the estimation of factor matrices from an input-output tensor
model. However, it embeds temporal dynamics in a higherorder tensor whereas the parametrization we propose follows
the control engineering approach to combine the temporal dynamics linearly while modeling independently each coefficient
matrix with a sum of Kronecker matrices. Besides, we allow
the Kronecker rank to be strictly larger than one for more
generality and applicability for identification and control of
systems such as adaptive optics. These two points are crucial
to achieve good accuracy estimations in e.g a laboratory
environment and hence, enable its effective use for control.
Third, the QUARKS methodology proposes regularization to
estimate stable and sparse models.
Another work related to the framework we propose deals
with blind source separation using tensor representations, [24].

The approach consists in estimating two matrices M and S
from the measurements stored in X given the relationship:
X = MS

(3)

where M represents the mixing matrix and S ∈ Rn×K
the n source signals for K time samples. The work [24]
relies on a low-rank decomposition of a certain reshaping
(equivalently, segmentation) for either/both the rows of the
mixing matrix and the source channels in order to achieve
a trade-off between data compression and accuracy of the
data fit. Both the present paper and [24] reshuffle the mixing
vectors/coefficient matrices in order to exhibit a low-rank
matrix and consequently, reduce the number of modeling
parameters. Nonetheless, our modeling assumptions differ in
three ways. We model the coefficient-matrices with lowerdimensional matrices without making restrictive assumptions
on the signals rather than being obtained from a regular grid
and being persistently exciting. We focus on the specific case
where the sources signals S are known which allows to get rid
of the ambiguity transformation inherent to BSS identification
and to formulate spatial and temporal stability constraints on
the coefficient-matrices Ai . Last, we exploit the 2D structure
of the network and separability of the modeled functions
in order to reduce the number of parameters. This point is
detailed in Section III. The different modeling assumptions
lead to distinct optimizations procedures.
In the following, the class of low-Kronecker rank matrices
is studied with a focus on modeling 2D spatial-temporal
dynamical systems of Vector Auto-Regressive form. The Kronecker tool as presented in this paper is meant to break down
the curse of dimensionality when working with arrays of
higher dimensions and without necessarily enforcing a priori
a sparsity pattern in the network, hence allowing to discover
both spatially varying dynamics and an unknown topology
from the data. It also serves as the basis for other more useful
identification approaches such as subspace identification, see
e.g [25]. As such, it will establish the fundamentals of a
new modeling framework for the identification and analysis
of large-scale 2D dynamical systems. The challenge lies in
deriving algorithms that are, on the one hand, scalable in terms
of data storage as well as in terms of computational complexity
in identifying and using these models, e.g in subsequent
control design, and on the other hand, that still ensures similar
prediction performances compared to the unstructured leastsquares estimates. The main contributions of this paper are
the definition of a new class of dynamical systems -of low
Kronecker rank-, the formulation of a regularized cost function
for identification and the formulation of an Alternating Least
Squares algorithm that is proven to converge to the global
minimum with O(N 3 Nt ) computational complexity where Nt
is the number of temporal samples.
The paper has the following outline. Section II describes
the class of sums-of-Kronecker matrices, while Section III
associates a VAR model associated with network data. In
Section IV we describe regularization methods to emphasize
the identification of stable models both in time and space. We
study in Section V the Alternating Least Squares algorithm
with a focus on the conditions to ensure global convergence.

3

Missing data is handled in Section VI, both for a random
position of the missing entries or for different sensor network
architectures than a rectangular grid. The methods are then
illustrated in Section VII on a random low-Kronecker rank
VARX model and a practical scenario dealing with open-loop
identification of the atmospheric turbulence for adaptive-optics
purposes.
Notations. Scalars are denoted by lower or uppercase letters
or symbols. Vectors are written as boldface lower-case letters
such as x. The boldface is used to make a distinction between
indexing a set of vectors, such as x1 , x2 , and referring to the
elements of a single vector x ∈ Rn , such as x1 , . . . , xn . The
null vector and the vector of ones is denoted by 0 and 1
respectively, where an index can be used to explicitly show
its size e.g 1n ∈ Rnp
. The Euclidean norm of a vector x is
x21 + . . . + x2n = hx, xi. The sum in
written as kxk2 =
absolute P
value for the elements in x ∈ Rn is denoted with
n
kxk1 = i=1 |xi |.
Matrices are represented by boldface uppercase letters such
has X. The element located at the i-th row and j-column of
the matrix X is written as xi,j , or x?,(i,j) when the matrix
is denoted with X? . The inverse and transpose are written
as X−1 and XT respectively. The notation BDiag(Xi , N )
forms a block-diagonal matrix with X1 to XN located on
the block-diagonal. For a block-diagonal matrix X, the ith block is denoted with X [i]. M ATLAB-like notations are
used to denote columns and rows of matrices, e.g X(:, i)
refers to the i-th column of X, X(i, :) the i-th row. The
vectorization operator applied on X is written with vec(X) =

T
x1,1 x2,1 . . . xm,n . The operation of reshaping a vector into a matrix is denoted with ivec, e.g ivec(vec(X)) = X.
The Kronecker product of two matrices X, Y is represented by
the symbol ⊗ such as X⊗Y. The Frobenius
for a matrix
qP norm
m Pn
2
X ∈ Rm×n is denoted with kXk2F =
i=1
j=1 xi,j .
The big-O notation is used for describing computational
complexities and indicates the asymptotic growth rate of the
computational cost for a given mathematical operation. E.g
an operation costing O(n) floating-point operations (flops)
finishes in at most c · n flops, for some constant c.
Other section-specific notations are introduced in the respective section.

II. P RELIMINARIES
In this section, we review some of the most important
properties related to a low-Kronecker rank decomposition
of matrices. Such a decomposition relies on the existence
of block-matrices of equal size and that allow for a reorganization of the entries into a low-rank reshuffled matrix.
Definition 1. [14] Let m1 , n1 , m2 , n2 ∈ R. Let X ∈
Rm1 m2 ×n1 n2 and Xi,j ∈ Rm2 ×n2 such that:


X1,1
 ..
X= .
Xm1 ,1

···
..
.
···


X1,n1
.. 
. 
Xm1 ,n1

then the re-shuffle operator R(X) ∈ Rm1 n1 ×m2 n2 is defined
as:


T 
vec
X
1,1




..


.


T 



 vec X
m1 ,1





T 



 vec X1,2


R(X) = 

.
..




T 



 vec Xm ,2
1




..




.
 
T 
vec Xm1 ,n1
There exists a permutation matrix P in the set
Rm1 n1 m2 n2 ×m1 n1 m2 n2 such that:
vec(R(X)) = Pvec(X)

(4)

Lemma 1. [14] Let X be defined as in Definition 1, and let
X = F ⊗ G, with F, G ∈ Rm1 ×n1 × Rm2 ×n2 . Then:
R(X) = vec(F)vec(G)T
The operation in Lemma 1 can also be reversed by the
definition of the inverse vec operator ivec(.).
Lemma 2. [14] Let X be defined as in Definition 1, and let
an SVD of R(X) be given as:
R(X) =

r
X

σ` u` v`T

`=1



and let ivec u`



 
= U` , ivec v` = V` , then:
X=

r
X

σ` U` ⊗ V`

`=1

The integer r is called the Kronecker rank of X with respect
to the chosen block partitioning of X as given in Definition
1. The Kronecker rank of X is denoted with K-rank(X).
From Lemma 2, looking for a low-Kronecker rank approximation of a matrix is equivalent to finding a low-rank
approximation of the reshuffled matrix. The reshuffling operator R as defined in Definition 1 that yields a reshuffled
matrix of minimal rank r is not unique: reshuffling the blockmatrices row-wise rather than column-wise would yield the
same Kronecker-rank for X.
Definition 2. [α-decomposable matrices, [12]]
Let us consider a network of subsystems such that the latter
belong to α different classes, themselves composed of Ni
N ×N
subsystems.
be a pattern matrix. Define
Pj Let P ∈ R
βj =
N
(with
β
=
0)
and I[a1 :a2 ] as an N × N
i
0
i=1
diagonal matrix which contains 1 in the diagonal entries of
indices from a1 to a2 (included) and 0 elsewhere, then an
α-decomposable matrix (for a given α) is a matrix of the
following kind:
α 

X
M=
I[βi−1 +1:βi ] ⊗ L(i) + I[βi−1 +1:βi ] P ⊗ N(i)
i=1

4

The matrices L(i) are the diagonal blocks of M that model
the local dynamics, while the influence from the neighborhood
is represented by the matrices N(i) , according to the structure
of P.
When a state-transition matrix of a state-space model belongs
to the class of α-decomposable matrices, the associated network has a known interconnection pattern whose adjacency
matrix is P while α represents the number of non-identical
subsystems in the network. The pattern matrix P is allowed
to be time-varying. For α = 1 (and β1 = N ), these matrices
are simply called decomposable matrices.
As a generalization of this class of structured matrices, we
define next the class of sums-of-Kronecker product matrices.
Definition 3. The class of sums-of-Kronecker product matrices
contains matrices of the following kind:
M=

r
X

Ma (i) ⊗ Mb (i)

i=1

with Ma (i) ∈ Rm1 ×n1 and Mb (i) ∈ Rm2 ×n2 . This class is
denoted with Kr .
The matrices Ma (i) , Mb (i) are called factor matrices.
With this class of sums-of-Kronecker matrices, it is not
necessary to have knowledge of a pattern matrix P as with decomposable matrices. Therefore, the topology of the network
need not to be known in advance. Moreover, the network may
be composed of heterogeneous subsystems without any further
specifications on the structure of the factor matrices. When
describing large-scale networks, this structure is advantageous
for its high compression capabilities. While m1 m2 n1 n2 entries are necessary to describe M in the unstructured case,
only r(m1 n1 + m2 n2 ) elements are required in the Kronecker
framework. We first review a few results on the Kronecker
ranks before studying the computational complexity of elementary operations when the matrix belongs to Kr .
Lemma 3. For two matrices A ∈ KrA , B ∈ KrB such that
A=

rA
X

A`,i ⊗ Ar,i ,

i=1

B=

rB
X

•
•

The product AB reads:
AB

=
=
=

Proof. • The result readily follows from expressing A + B
in Kronecker format with:
A+B =
=

rA
X

A`,i ⊗ Ar,i +

i=1
rA
+rB
X

rB
X

A`,i ⊗ Ar,i

i=1
rA X
rB
X

rB
X

B`,j ⊗ Br,j

j=1

A`,i B`,j ⊗ Ar,i Br,j

i=1 j=1
rX
A rB

M`,i ⊗ Mr,i

(5)

•

for some matrices M`,i , Mr,i . There is a maximum of
rA rB terms in the sum.
2
We can now infer K-rank(A2 ) ≤ rA
, and it follows
i
i
K-rank(A ) ≤ rA for some integer i.

No theoretical result is however known (to the best of the
authors’ knowledge) on the Kronecker rank of A−1 when
r > 1. Lemma 4 provides insight on the benefits to use the
class of Kronecker matrices to speed up simple linear algebra
operations.
2
Lemma 4. Let x ∈ RN . Then, the orders of magnitude of the
computational complexity orders for matrix-vector multiplication, matrix-matrix multiplication and matrix inversion is as
follows:
2

Ax
AB
A−1 (case r = 1)

A, B ∈ RN ×N
O(N 4 )
O(N 6 )
O(N 6 )

2

A, B ∈ Kr
O(rN 3 )
O(r2 N 3 )
O(N 3 )

The complexity obtained with the Kronecker parametrization considers the operations required for forming the factor
matrices only.
Proof. From (5), the cost for computing the left factor matrix
Ml,i is O(N 3 ) for each 2r2 matrix-matrix products. Computing the Kronecker product M`,i ⊗Mr,i costs O(N 4 ) however.
The inverse for A is determined via:

B`,i ⊗ Br,i

K-rank(A + B) ≤ rA + rB
K-rank(AB) ≤ rA rB
i
For all i ∈ N, K-rank(Ai ) ≤ rA

rA
X

i=1

−1
A−1 = A−1
`,1 ⊗ Ar,1

i=1

where A`,i , Ar,i , B`,i , Br,i ∈ RN ×N , the following properties
hold:
•

•

−1
3
Computing A−1
`,1 and Ar,1 costs O(N ).

Approximating the inverse of large-scale low-Kronecker
rank matrices A ∈ Kr , r > 1 is an on-going research topic
which the works in [20] and [21] have investigated.
From Lemma 4, significant data compression is achieved when
r  N , which is the class of Kronecker models we are
interested in and is the focus of the following derivations.
III. P ROBLEM FORMULATION

B`,i ⊗ Br,i

i=1

Low-Kronecker rank matrices are now used to model the
input-output relationship of 2D networked systems.

C`,i ⊗ Cr,i

i=1

where C`,i = A`,i if i ≤ rA and B`,i otherwise, and
similarly for Cr,i with the right factor matrices.

A. QUARKS models
Let us consider a network with N ×N nodes, each of which
is associated with a scalar sensor signal. In the sequel, we

5

assume N > 1. The sensor readings at the time instance k are
stored in the matrix Sk ∈ RN ×N as:


sk,(1,1) sk,(1,2) · · · sk,(1,N )
 sk,(2,1) sk,(2,2)
sk,(2,N ) 


(6)
Sk =  .

.
..
..
..
 ..

.
.
sk,(N,1)

sk,(N,2)

···

sk,(N,N )

Although the framework that we present here extends straightforwardly to arrays with nodes having multiple outputs, we
only dwell on this case in Section VII. In this paper we will
consider that the temporal dynamics of this array of sensors
is governed by the following VAR(X) model:
p
  X


 
vec Sk =
Ai vec Sk−i + vec Ek
(7)
i=1





with vec Ek zero-mean white noise with identity covariance matrix. Covariance estimation for low-Kronecker rank
matrices has been addressed in [17] and is not the subject of
further investigations in this paper. The coefficient matrices
Ai in the VAR model (7) are in general highly structured. We
consider the case they belong to the set Kr and focus on the
coefficient matrices Ai . To address an identification problem
we parametrize these coefficient matrices as:
Ai =

ri
X

(j)

(j)

Ai ,

Ai

(j)

(j)

= M(bi )T ⊗ M(ai )

(8)

j=1
(j)

(j)

with the vectors ai and bi parametrizing the matrices
(j)
(j)
M(a
 i ) and M(bi ) in an affine manner. With the notation
vec Sk = sk , the VAR model (7) can be rewritten as,

B. The identification problem of QUARKS models.
Given the model structure of QUARKS models, the problem
of identifying these models from measurement sequences
t
{Sk }N
k=1 is fourfold:
1) The temporal order index p.
2) The spatial order index ri for each coefficient matrix.
(j)
3) The parametrization of the matrices M(ai ) and
(j)
M(bi ). An example of a parametrization of the ma(j)
(j)
trices M(ai ) and M(bi ) is (block) Toeplitz.
(j)
(j)
4) The estimation of the parameter vectors ai , bi up to
an ambiguity transformation. This requires the specification of a cost function. An example of such a cost
function using the model (10) is the following least
squares cost function,
p
Nt


X
X
kSk −
Mai Iri ⊗ Sk−i Mbi k2F (12)
min
(j)

(j)

ai ,bi

k=p+1

for data batches with Nt points.
By the selection of the parameter p and the particular
choices of the parametrization in step 3 above, various special
cases of restricting the coefficient matrices Ai in (7) to particular sets such as Kri can be considered. Further constraints
to the least-squares cost function (12) might be introduced to
(j)
(j)
look for sparsity in the parametrization vectors ai and bi .
The non-uniqueness of the optimal solution for the cost
function (12) is highlighted next. One way to solve this
estimation problem is via vectorization of the sensor signals
Sk :

sk

=

i=1

(j)
M(bi )T

⊗

(j)
M(ai )


sk−i + ek (9)

Ai sk−i k22

i=1
ri
X

(j)

(j)

M(bi )T ⊗ M(ai )

From (4), the reshuffling operator R(.) is a bijection in
2
2
RN ×N , therefore the above minimization problem is equivalent to:
min
(j)

(10)

s.t



(1)

···


(r )
M(ai i )

k=p+1

ksk −

p
X

Ai sk−i k22

i=1

R(Ai ) = Ui ViT

(13)

where:
Ui
(11)

where
(1)
M(ai )

Nt
X

(j)

ai ,bi

i=1

=

p
X

j=1

This can also be written explicitly as,
p


X
Sk =
Mai Iri ⊗ Sk−i Mbi + Ek

Mbi

Ai =

j=1

h

ksk −

k=p+1

s.t

j=1

we can write the VAR model (9) as,
p X
ri

X
(j)
(j)
Sk =
M(ai )Sk−i M(bi ) + Ek

Mai

(j)

ai ,bi

Using the following Kronecker rule, for matrices X, Y, Z of
compatible dimensions such that the product XYZ exists,

  


ZT ⊗ X vec Y = vec XYZ

i=1

Nt
X

min
(j)

p X
ri
X

i=1

i

M(bi )


..

= 
.


(ri )
M(bi )

The VAR(X) models (9), (10) or (11) are called Kronecker
VARX network models and abbreviated with QUARKS models.

Vi




i
(1)
(r)
vec M(ai )
. . . vec M(ai )
h 


i
(1)
(r)
= vec M(bi )
. . . vec M(bi )

=

h

For a non-singular transformation Ti ∈ Rr×r , the constraint
(13) can be equivalently written as:
R(Ai ) = Ũi ṼiT
T
where: Ũi = Ui Ti , ṼiT = T−1
i Vi . The non-uniqueness of
the factor matrices is not an issue for practical use of QUARKS
models as it does not affect the prediction-error.
Remark 1. Let m ∈ {1, .., N 2 }. Blind source separation (3) as
described in [24] reshapes either (or both) the mixing vectors

6

M(m, :) and sources S(m, :) in (3) to form low-rank matrices.
Then, there exists different left and right matrices for each
T
mixing vector M(m, :) such that R(M(m, :)) = um vm
, or
equivalently,
M(m, :) =

r
X

um (:, j)T ⊗ vm (:, j)T

j=1
I×r

where um ∈ R , vm ∈ RJ×r for two scalars I, J. The
parameters I, J are user-defined contrary to the QUARKS
modeling, where I, J = N . Hence, all mixing vectors are
decoupled independently contrary to the description for the
QUARKS model (8) that considers a low-rank segmentation
of both the rows and columns of the mixing matrix M.
We illustrate in the case where p = 1 and M = A1 .
If rank(R(M)) = r, then rank(R(M(m, :))) = r and
rank(R(M(:, m))) = r. Fixing I, J to N and considering N 2
sources, there are 2rN 3 unknown coefficients to estimate while
the modeling (8) represents the coefficient matrices with 2rN 2
entries. The QUARKS modeling decreases the data storage
requirements by an order of magnitude.
An important challenge in solving the parameter estimation
problem (12) is the computational efficiency for the case when
the size N of the array is assumed to be large. A solution to
(13) based on nuclear norm minimization was proposed in
[26] but does not allow for an efficient integration of stability
constraints.
IV. R EGULARIZATION INDUCING SPATIAL - TEMPORAL
STABILITY AND SPARSITY

The Kronecker rank is assumed equal for all i, i.e ri = r,
without constraining the insights in this section.
A. Stability of VAR models
A first solution to enforce stability for VAR models of the
form (7) is to formulate the Lyapunov equation associated with
the model in companion form. Denote the companion matrix:


A1 . . .
. . . Ap
IN 2 0
0



.. 
.
.
..
..

.
Ā = 




.
.
..
..


IN 2 0

Nonetheless, it is not suitable for large VAR models as the
dimension of the LMI scales with 2N 2 , which is prohibitive
for large-scale networks. The Kronecker structure does not
help the formulation of an efficient algorithm for two reasons.
First, the properties of the Kronecker product [14] are onesided, e.g for two matrices A, B, if A and B are positivedefinite, then A ⊗ B is positive-definite while the reverse is
not true in general. Second, there is no relevant literature for
efficiently solving large-scale LMI with low-Kronecker rank
matrices.
In [7], the stability for VAR models is guaranteed by modeling
the impulse response from one node to all the other ones
in the network as a zero-mean Gaussian process and with
an adequately chosen covariance matrix, which ensures that
the parameters of the impulse response are decaying with
increasing temporal index. We refer to [27] for a general
introduction to kernel methods and to [28] [29] for an application to system identification. In the following paragraph
we integrate these results as an additional regularization to the
cost function (12). We introduce the positive-definite matrix
Pt ∈ Rp×p following a Gaussian-kernel to fit stable impulses.
For example, a Diagonal-Correlated (DC) kernel Pt is defined
with:
i+j
pt,(i,j) = ξ 2 η |i−j|
for i, j = 1..p, and where the optimal hyperparameters
−1 ≤ η ≤ 1, 0 ≤ ξ < 1 shall be determined either by grid
search or within the framework of Bayesian optimization to
tune both the decay rate and the smoothness of the impulse
response. Let Wt be a square root of P−1
t . As there is no prior
information nor physical meaning to distinguish between the
different factor matrices, these are regularized independently
with the cost:


U1 (:, j)V1 (:, j)T
r
X

 2
..
kQt 
 kF
.
T
j=1
Up (:, j)Vp (:, j)
where Qt = Wt ⊗IN 2 . In a more compact notation, we write:


U1 (:, j)V1 (:, j)T


..
f (Ma (j) , Mb (j) ) = 

.
rt (Ma , Mb )

The matrix inequality originating from Lyapunov theory reads:
T

P − ĀPĀ  0,

P0

where P is symmetric positive-definite. Rewriting this Linear
Matrix Inequality (LMI) using the Schur complement yields:


P Ā
 0,
P  0,
PQ = I
(14)
ĀT Q
The following cost-function does ensure the identification of
a stable VAR model,
min

Ai ,P,Q

s.t

Nt
X
k=p+1

(14)

ksk −

p
X
i=1

Ai sk−i k22

Up (:, j)Vp (:, j)T
r
X
=
kQt f (Ma (j) , Mb (j) )k2F
j=1

Such a regularization rt (.) is bilinear in the unknowns
Ma , Mb .
B. Spatial sparsity
Real graphs or the regular networks from discretized Partial
Differential Equations are such that each node is connected
to a very limited number of other nodes with respect to
the network’s size. In the latter case, the neighborhood is
localized which gives rise to a multi-banded structure of the
full coefficient-matrices, equivalent to a banded structure of
each factor matrix. In case of high coupling, as is observed e.g
in the atmospheric turbulence modeling as discussed in Section

7

VII, we rather tune the decay of the parameters away from the
main diagonal rather than minimizing the number of non-zero
entries. Furthermore, it will become clear in the next section
that we would like to avoid all non-differentiable functions in
the cost function, hence the focus is laid on kernel methods
rather than on minimizing the `1 -norm of the factor-matrices.
In this framework, an exponentially decreasing sequence has
been studied in [7] for sparse network identification.
For a given matrix X ∈ RN ×N , denote (for i ≥ 0):


diag(X, i)T = x1,1+i x2,2+i . . .


diag(X, −i)T = x1+i,1 x2+i,2 . . .
d0
∀i = 1..N − 1, di

= diag(X, 0)

= diag(X, i)T

diag(X, −i)T

with Xi ∈ RN ×N . Reshuffling the matrix X according
to Definition 1 is rank(R(X)) ≤ N .
2) Consider e.g




0 1
4 5
X=
,
Y=
2 3
6 7
The matrix X ⊗ Y is not block-Toeplitz neither has
Toeplitz blocks but has Kronecker rank 1.

T

We reshape the elements of a square matrix diagonal-wise,
starting by the main diagonal, and denote this operation with
the operator D:

T
2
D(X) = dT0 dT1 . . . dTN −1 ∈ RN
2

1) Let X be symmetric block-Toeplitz such that:


X1 . . . XN

.. 
..
X =  ...
.
. 
XN . . . X1

Proof.

2

Let us introduce a diagonal matrix Ks ∈ RN ×N such that:


IN k 1
0
...
0

.. 
..
 0
.
I2(N −1) k2
. 


Ks =  .

.
.
.
.
.
 .
.
.
0 
0
...
0 I2 k N
where the scalars ki are such that 0 < ki < ki+1 . For example,
a valid choice of such scalars is ki = eζi with ζ > 0. The prior
(j)
(j)
for matrices M(ai ), M(bi ) with values decaying away
from the main diagonal is then:
p X
r

 
T
X
(j)
(j)
kKs D M(ai ) D M(bi ) KTs k2F
i=1 j=1

This spatial regularization is denoted with rs (Ma , Mb ).
C. Structured factor matrices
The parametrization of the factor matrices based on prior
knowledge of the network may help either to further reduce
the computational complexity of the model identification step,
or to cast the model into a structure useful for control. The first
category include banded, symmetric, Toeplitz and circulant
patterns. Exploring such structures on the factor matrices is
very attractive numerically as the number of parameters to be
estimated reduces further. The block-Toeplitz Toeplitz-blocks
(bTTb) structure arises e.g when modeling 2D homogeneous
spatially-invariant phenomena on a rectangular grid. Many
functions in optics are isotropic, for example the Point Spread
Function or covariance matrix of the atmospheric turbulence,
and can be modeled with a sum of few Kronecker terms. The
Kronecker and bTTb structures are related, but not equivalent.
2
2
Lemma 5. Let X ∈ RN ×N . The following holds:
• If X is symmetric block-Toeplitz, then X has a Kronecker
rank of at most N .
• If X has a Kronecker rank of 1, it does not in general
imply neither that X is block-Toeplitz nor has Toeplitzblocks.

The second category contains for example the sparse (with
unknown pattern) or SSS structure. The SSS structure is
different from the Toeplitz, and has advantages especially
when it comes to model spatially-varying systems. As mentioned in Section I, the efficient use of SSS matrices has
been thoroughly studied in [8] while the extension to MultiLevel structures is an on-going research question. Modeling
each factor matrix of the model as SSS enables significant
improvements in the computational cost for future simple
linear algebra operations. For example, the cost for standard
matrix computations (×, .−1 ) scales linearly with respect to
2
2
the matrix size. Inverting a matrix M ∈ RN ×N written as
M = M1 ⊗ M2 in which both M1 , M2 have a SSS structure
requires O(N ) operations instead of O(N 6 ). Because such
(j)
(j)
a parametrization for the matrices M(ai ), M(bi ) is not
(j)
(j)
affine in the parameters ai , bi , the identification of the
SSS matrices is performed offline, i.e after having obtained
(j)
(j)
an estimate for the non-parametrized M(ai ), M(bi ).
D. The regularized cost function for QUARKS identification
The cost function for the identification of sparse stable
QUARKS models reads:
p
Nt


X
X
min
kSk −
Mai Ir ⊗ Sk−i Mbi k2F
(j)

(j)

ai ,bi

k=p+1

i=1

+µrt (Ma , Mb ) + λrs (Ma , Mb )

(15)

where µ, λ are regularization parameters. The cost function
(15) belongs to the class of multi-convex problems in which
fixing one set of variables yields a convex problem. Adding
regularization to the cost function aims at decreasing the
prediction error of the estimated VAR model when dealing
with noisy and short data batches rather than speeding up the
convergence as done in [35].
Remark 2. The regularization in (15) is bilinear contrary
to the one analyzed in [30], [31] within the framework
of Principal Component Analysis (PCA). Based on [30], a
regularization for the minimization (13) could be proposed as
follows:
p
Nt


X
X
min
ksk −
Ai sk−i k22 + µ kUi k2F + kVi k2F
Ai ,Ui ,Vi

s.t

k=p+1

i=1

∀i = 1..p,

R(Ai ) = Ui ViT

(16)

8

The regularization does not independently minimize kUi k2F or
kVi k2F while it handles directly the elements in the coefficient(j)
matrices Ai in (15).
V. B I - CONVEX C OST FUNCTION APPROACH
A. An Alternating Least Squares approach
The regularized least-squares representation (15) is bilinear
in its unknowns but features relatively small factor matrices,
which has the advantage that constraints on the parametrization
(j)
(j)
of the matrices M(ai ) and M(bi ) can be more easily taken
into consideration than via a low-rank minimization on the
large-scale reshuffled matrix, [26]. A non-linear optimization
scheme such as the separable least-squares [34] proceeds
with two steps, one of which however consists of non-linear
optimization. Iterative algorithms have been derived as a
generalization of the linear Gauss-Seidel iterations for solving
coupled Sylvester matrix equations in [33]. Similarly as in
[23], we propose to address (13) by solving a sequence of
linear least-squares and using ALS, which is a special case of
the block non-linear Gauss-Seidel method as highlighted in
[35].
(j)
(j)
We denote the κ-th column of M(ai ) and M(bi ) with
(j)
(j)
respectively ai,κ and bi,κ . The data-fitting term in (15) is first
rewritten with:
kS̃ − M̄a Xb k2F

where:
√
Gb (Ma )
Gb,j (Ma )

=

λPr,s BDiag(Gb,j (Ma ), r)Pc,s
(j)

= BDiag(Ks ⊗ Ks D(M(ai ), p)

The matrices Pr,s and Pc,s permute respectively the rows
and columns such that Gb (Ma ) is block-diagonal. The cost
function is then separable for each column of Xb :
Gb (Mb (`) )xb =

N
X

Gb (Mb ) [i] Xb (:, i)

i=1

Therefore, the cost function (15) for updating Xb reads now:
  

N
M̄a
S̃
X
min
k  0  −  Fb (Ma )  Xb (:,i)k2F
(17)
Xb
i=1
Gb (Ma ) [i]
0
At iteration `, the solution is denoted with Xb (`) . Similarly,
the least-squares for updating Xa is:
  

N
M̄b
S̃
X
min
k  0  −  Fa (Mb )  Xa (:,i)k2F
(18)
Xa
i=1
Ga (Mb ) [i]
0
where:

where:


(j)
IN ⊗ bi,1


..

(IN ⊗ Ũi ) 
.


(j)
IN ⊗ bi,N


M̄b,1,1 . . . M̄b,1,r . . . M̄b,p,r


Ma1
 .. 
xa = vec(Xa )
 . ,


S̃j,i




sp+1,(j,i)
S̃1,1



.
.
..
=
 , S̃ =  ..


S̃1,N
.. 
. 

...

sNt ,(j,i)
S̃N,1 . . . S̃N,N


Sp+1−i (1, :) . . . Sp+1−i (N, :)


..
..
Ũi = 

.
.
SNt −i (1, :)

...


M̄a,i,j

=

M̄a,i

=

M̄a

=

Xb

=

SNt −i (N, :)
(j)
ai,1

M̄b,i,j

=

M̄b

=

Xa

=

Map



⊗ IN


..


(IN ⊗ Ũi ) 
.

(j)
ai,N ⊗ IN


M̄a,i,1 . . . M̄a,i,r


M̄a,1 . . . M̄a,p


Mb1
 .. 
xb = vec(Xb )
 . ,
Mbp

The regularization µrt (Ma , Mb ) is also written as
kFb (Ma )Xb k22 , where each block Fb (Ma ) [i, j] concatenates
in the block-diagonal the terms ranging from 1 to r:


Fb (Ma ) [1, 1] . . . Fb (Ma ) [1, p]


..
..
..
Fb (Ma ) = 

.
.
.
Fb (Ma ) [p, 1] . . . Fb (Ma ) [p, p]
√
(j)
µwt,(i,j) BDiag(IN ⊗ vec(M(ai )), r)
Fb (Ma ) [i, j] =
Moreover, a matrix Gb (Ma ) is derived such that
λrs (Ma , Mb ) = kGb (Ma )xb k22

The solution is denoted with Xa (`) .
The least-squares (17) and (18) are iteratively solved starting
(0)
with some random initial guess for Ma until some stopping
criterion is reached. Denote the value of the cost function
(18) with c(`) . The iterations are stopped when the relative
decrease of the cost function is lower than a given threshold
: |c(`) − c(`−1) | ≤  during Imin iterations. Algorithm 1 summarizes the steps. Normalization of the columns Xb (`) (:, i)
for i in the set {1, .., N } along with the scaling respectively
in the lines 10-12 and 28-31 are studied in Algorithm 1. It
is considered within the scope of this paper for two reasons.
First, this normalization is a key ingredient in deriving the
convergence proof to the global minimum in V.C. More
precisely, it guarantees the existence of a unique fixed point.
Second, the non-uniqueness of the solution may imply that
the estimated factor matrices scale to very large (resp. very
low) values, which the normalization prevents from happening.
It however requires the knowledge of Xb (:, i) which may
not be always available in practice. Numerical examples in
Section VII.A illustrate the impact of normalization on the
convergence.

9

Algorithm 1 ALS for QUARKS identification
1: procedure ALS({Sk }, r, p, {`max , , Imin }, {η, ξ, ζ, λ, µ},
{kXb (:, i)k2 }1..N )
(0)
2:
Initialize
Ma = randn(Nrp,N), ` = 0, I = 0.
 
S̃
3:
Y :=  0 .
0
4:
while ` < "`max and I < #
Imin do
(`−1)
M̄a
5:
F0 :=
.
Fb (Ma (`−1) )
6:
Compute FT0 F0 .
7:
for i = 1..N
 do

F0
8:
F :=
.
Gb (Ma (`−1) ) [i]
9:
Xb (`) (:, i) := (FT F)−1 FT Y(:, i)
10:
if isempty({kXb (:, i)k2 }1..N ) = 0 then
kXb (:,i)k2
11:
Xb,n (`) (:, i) = Xb (`) (:, i) kX
(`) (:,i)k
2
b
12:
end if
13:
end for"
#
(`)
M̄b,n
14:
F0 :=
.
Fa (Mb,n (`) )
15:
Compute FT0 F0 .
16:
for i = 1..N
 do

F0
17:
F :=
.
Ga (Mb,n (`) ) [i]
(`)
18:
Xa (:, i) := (FT F)−1 FT Y(:, i)
19:
end for
(`)
20:
c(l) := kFXa − Yk2F .
(`)
(`−1)
21:
if |c − c
| <  then
22:
I := I + 1
23:
else
24:
I := 0
25:
end if
26:
` := ` + 1
27:
end while
28:
for i = 1..p do
(`)
(`)
d
29:
M
bi = Mb i sign(mb,i,(1,1) )
(`)
(`)
d
30:
M
ai = Ma i sign(mb,i,(1,1) )
31:
end for
32: end procedure

B. Computational complexity
For matters of comparison, we evaluate first the complexity
for estimating the coefficient-matrices associated with the
unstructured VAR. Using (7) with temporal data within the
range {1, ..., Nt } with Nt ≥ N 2 p to recover a unique solution,
we write:


s̃f = A1 . . . Ap s̃p + ẽp

Operation
s̃p s̃Tp
(s̃p s̃Tp )−1
s̃f s̃Tp

Flops
O(N 4 Nt )
O(N 6 )
O(N 4 Nt )

The dependency on the number of temporal samples is kept:
a correct identification in noisy conditions often requires
Nt  N 2 p. Hence, the complexity for estimating unstructured
VAR is O(N 4 Nt ) where Nt  N 2 p.
The lines 6, 9, 15 and 18 are the most computationally costly
of Algorithm 1. There are two case worth investigating: λ = 0
and λ" 6= 0. First, if #λ = 0, the pseudo-inverse for the matrix
(`−1)
M̄a
F=
is computed only once at each iteration.
Fb (Ma (`−1) )
(`−1)

requires (Nt − p)rp matrix-matrix multipliForming M̄a
cations of size N × N . The number of temporal samples is
such that N (Nt − p) ≥ N rp to guarantee a unique solution
of each subproblem without regularization. In practice, larger
sizes of network require a larger number of temporal samples,
and therefore Nt is included in the computational count as
it depends on N . Assuming that r is independent of N , the
complexity is O(N 3 Nt ) flops. Moreover, only X(Ma (`−1) ) is
dense, therefore computing FT F costs O(N 3 Nt ). Its inverse
requires O(N 3 ), while right-multiplying the latter with FT
reaches O(N 3 Nt ). Assuming the number of ALS iterations
does not depend on N , the computational complexity for
Algorithm 1 with λ = 0 reaches O(N 3 Nt ) where Nt  rp.
If λ 6= 0, the matrix F is partitioned in 3 parts. Computing
FT F in line 9 is negligible because Gb (Ma (`−1) ) [i] is sparse.
Moreover, the term FT0 F0 is computed once in line 7 with
O(N 3 Nt ) flops. Computing the inverse of FT F requires
O(N 3 ) flops, while multiplying the inverted matrix with FT
costs O(N 3 Nt ). These two operations need to be repeated N
times, although it should be performed in parallel. The price
for computing the lines 15 and 18 is similar to the above
discussion. When the algorithm is computed sequentially and
without making use of distributed computing platforms, the
overall cost reaches O(N 4 Nt ).
The costs are summarized below.
Operation
Lines 6 and 15
Lines 9 and 18 (for each i)
Lines 9 and 18 (total for all i)
Total (with λ 6= 0)
Total (with λ = 0)

Flops
O(N 3 Nt )
O(N 3 Nt )
O(N 4 Nt )
O(N 4 Nt )
O(N 3 Nt )

C. Convergence proof for the normalized Alternating Least
Squares

where


s̃f = sp+1

The least-squares estimation for the coefficient-matrices is
hence equal to:
h
i

−1
c1 . . . A
cp = s̃f s̃Tp s̃p s̃Tp
A

...


sNt ,


sp
 ..
s̃p =  .

...


sNt −1
.. 
. 

s1

...

sNt −p

In this paragraph, we consider the normalized version of
Algorithm 1 and prove its convergence to a global minimum.
We assume the temporal order p and spatial order r to be
both equal to 1 and λ = 0. Further remarks on more general

10

spatial orders are made at the end of Section V. As noted in
the previous paragraph, the algorithm simplifies when λ = 0.
The updates for minimizing the regularized cost function (15)
now read:
  

M̄a
S̃
Xb k2F
(19)
min k
−
Fb (Ma )
0
Xb
|{z} | {z }
¯
M̄
a

Ỹ




M̄b
kỸ −
Xa k2F
Fa (Mb )
| {z }

min
Xa

(20)

¯
M̄
b

We denote with a the vector concatenating all columns ai
of Ma (:= Xa ). The estimate of a at iteration ` is denoted with a(`) . Similar notations hold for Mb (:= Xb when
(p, r) = (1, 1)). The convergence proof relies on the work
in [36] that establishes the result when Ma , Mb are vectors
and without regularization. Therefore, we review the results
in the following for completeness and highlight the nontrivial extensions in Appendix A. The convergence proof
uses the Contraction Mapping Theorem, [37]. A functional
representation of the three steps in Algorithm 1 reads:
b(`)
bn (`)
a(`)

= F1 (a(`−1) )

a

Proof. The proof is derived in Appendix A.
The assumption A2 corresponds to the persistency of excitation from the data and is a key ingredient in the convergence. Recasting the ALS into an input-output representation
and assuming A2 has enabled to avoid rank-deficiencies in
the matrix F and therefore swamps as observed for tensor
decomposition in [35] do not occur.
Theorem 2. Let the fixed point from Theorem 1 and assume
cb to the least-squares
A2. Then the solution X


M̄a
min kỸ −
Xb k2F
Fb (Ma )
Xb
is unique.
¯ is full column
Proof. The solution is unique if the matrix M̄
a
rank. The Gramian is written as:
¯ T M̄
¯ = M̄T M̄ + F (M )T F (M )
M̄
a
a
a
b
a
b
a
a
The first term is detailed into:

= F2 (b(`) )
= F3 (bn (`) )

M̄Ta M̄a =

These equations can be expressed using a single operator F(.)
mapping the estimate a(`−1) to a(`) :
(`)

Then, the map F : Xa → Xa is a contraction on Xa when
b ∈ Xa .
Nt → ∞ and has a unique fixed point a

= F3 (F2 (F1 (a
(`−1)

= F(a

(`−1)

)))

p X
r X
N
X

(j) T

(j)

(ai,n ⊗ IN )ŨT Ũ(ai,n ⊗ IN )

i=1 j=1 n=1

With the assumption A2, (ŨT Ũ+µIN 2 ) is full rank N 2 , and
¯ T M̄
¯ is full rank, N .
when Ma is non-zero, M̄
a
a

)

Lemma 6. [The Contraction Mapping Theorem, [37]] Let
(X, D) be a non-empty complete metric space where D is
a metric on X. Let F : X → X be a contraction mapping
on X, i.e., there is a nonnegative real number Q < 1 such
that D(F(x), F(y)) ≤ QD(x, y), for all x, y ∈ X. Then the
map F admits one and only one fixed point x? ∈ X which
means x? − F(x? ) = 0. Furthermore, this fixed point can be
found from the convergence
of an iterative sequence defined

by x(`+1) = F x(`) for k = 1, 2, ... with an arbitrary starting
point x(0) in X.
We define now two sets, associated to each of the variables b
and a:
Xb

2
= {b
b ∈ RN |∀i ∈ {1, ..., N }, kbbi k2 = kbi k2 , bb1 > 0}

Xa

= {b
a ∈ RN |∀i ∈ {1, ..., N }, kabi k2 ≤ kai k2 }

2

Theorem 1. Let p = 1, r = 1.
If the following statements are true:
• A1: the noise components in vec(Ek ) are independent
identically distributed (i.i.d) with zero-mean and finite
variance.
• A2: the matrix Ũ is full column rank.
• A3: either kbi k2 or kai k2 is known for all i and the first
non-zero entry of b or a is strictly positive.
(0)
• A4: the initial estimate Ma
is non-zero.

The Theorems 1 and 2 prove that, whatever the (nonzero) initial conditions are, the Alternating Least Squares in
Algorithm 1 converges to a global minimum asymptotically
when Nt → ∞. When the temporal order is strictly larger
than 1, the solution to an update in line 9 or 18 in Algorithm
1 is unique if and only if the matrix F is full column rank.
This condition provides with indications on how to choose
the initial guess. In practice, we choose randomly generated
initial guesses independent for each factor matrix such that
X(Ma (0) ) is full column rank with probability 1.
Remark 3. Higher Kronecker ranks. For both (p, r) strictly
greater than one, global convergence has been observed for
all numerical experiments carried out with the non-normalized
version of Algorithm 1, i.e skipping lines 10 to 12, if the
assumptions A1 and A2 are satisfied and if the initial guess
Ma (0) is full rank. If r > 1, we have observed numerically
that the normalized version of Algorithm 1 requires a lot more
iterations to converge to the global minimum than without
normalization.
Remark 4. Initialization. We highlight the importance of
choosing random initial guesses to converge to the global
minimum. For illustration we choose p = 1 and r > 1.
(0)
Then, for α ∈ N, α < r, if the rank of M̄a is αN ,
we have observed that the solution to Algorithm 1 without
normalization corresponds to the same solution that would
have been obtained by choosing rather r = α.

11

VI. G ENERALIZATIONS

10 5

Practical scenarios also exists for which the QUARKS
identification do not apply as presented in the previous sections
due to a non-rectangular architecture of the sensor nodes. In
this section we discuss on how to deal with other relevant
sensor architectures.

10 0

10 -5

10 -10

A. Non-rectangular measurement grids
In this paragraph, we aim at generalizing the Kronecker networks to non-rectangular measurement grids. For example, the
use of control methods for high resolution imaging requires the
modeling of 2D spatial-temporal dynamical systems. When the
sensor measurements are obtained from a rectangular array, the
coefficient-matrices of the VAR model are composed of blocks
of equal size such that a sums-of-Kronecker decomposition is
possible.
When only sensor measurements within a circular aperture
are available as seen in Figure 1 (left), liting the sensor
measurements only located within the aperture results in 2level coefficient-matrices in which the blocks have different
sizes. Recovering a QUARKS model from data sampled in
a circular array requires further information on the network
such as sparsity. When considering the rectangular embedding
of the aperture, the sensor measurements lying outside the
circular aperture can be set to 0, which has a negligible impact
on the overall quality of the QUARKS estimation if each node
has a localized influence on its neighborhood.
When the structure is hexagonal, numbering the nodes for
lifting the data such that the matrices Ai in (11) have a lowKronecker rank structure can be done according to the patterns
in Figure 1 (right). It guarantees that Ai is 2-level matrix
with blocks of equal size. The low-Kronecker rank structure
is present in the coefficient-matrices if the dynamical process
modeled is separable (or approximated as) in the new set of
directions, that is, {red, green} in Figure 1 (right).

Figure 1. Left: The data is available within the purple array, and the
QUARKS models are computed over the whole square aperture. Middle:
Some numbering patterns when the nodes in the network are organized with
a hexagonal structure.

10 -15
0

50

100

150

200

250

300

350

400

Figure 2. Singular values of the reshuffled matrix R(M) when considering a
noise-free input-output model y = Mu. Two 2D arrays are created: the first
one represents the coordinates (ai , bi ) used as inputs for a function f , the
second one represents the sampling points (xi , yp
i ) of f . The function f is
defined with fai ,bi (xi , yi ) = exp(−d/σ 2 ), d = (xi − ai )2 + (yi − bi ).
For each (xi , yi ), a N × N matrix is built which is lifted to form M(i, :).
Rect2Rect means the inputs and outputs grids are identical and with sampling
points regularly spaced, Hex2Rect means that the input grid is hexagonal while
the output grid is rectangular. Random2Random has the topology displayed
on Figure 1 (Right) for both inputs and outputs grids.

the inputs- and outputs- grid on which it is sampled. Another
type of irregularity in the sensor matrix Sk with its measurements sampled on a regular rectangular grid is when there
are missing sensor channels. Consequently, the coefficientmatrices in the lifted VAR model (ignoring the missing entries)
do not retain the Kronecker structure. The approach proposed
here is to lift the data by keeping the missing entries but
consider them as unknown. In this paragraph, we study the
estimation of the factor matrices as well as the missing entries
by working on the rectangular embedding when the position
of the failing sensors does not change over time. Estimation
of such models is very specific to the QUARKS framework
in so far as these entries are either left out with unstructured
VARX estimation, or mingled with blind source separation
[24] to deal with the bilinearities between the coefficients and
the data.
Denote the set of indices corresponding to known, resp.
unknown, entries, in Sk with Ω, resp. Ω̄. When the matrices
Ai are known, it suffices to formulate a linear least squares to
retrieve the missing data for all time samples. It is formulated
in the next lines. The VARX model is decomposed between
known and unknown data:
p
X
sk,Ω + sk,Ω̄ =
Ai (:, Ω)sk−i (Ω) + Ai (:, Ω̄)sk−i (Ω̄) (21)
i=1
2

Figure 2 provides with some intuition on how to choose the
spatial order r. The Kronecker rank r is related to the spatial
separability of the 2D discretized process. Any separable
functions in both horizontal and vertical coordinates, isotropic
functions such as the Bessel or Gaussian function as illustrated
in Figure 2, are meaningful candidates for a low-Kronecker
rank decomposition.
B. Missing data
We have seen in the previous paragraph that an efficient
QUARKS modeling depends on both the sampled function and

where sk,Ω ∈ RN contains 0 at the position of unknown
entries, and sk,Ω̄ = sk − sk,Ω . It is rewritten into a leastsquares cost function that estimates the missing data:
min
xk

Nt 
X

kdk − Fxk k22 + βkxk k22



k=p+1

where β is a regularization parameter and
dk

= sk,Ω −

p
X

Ai (:, Ω)sk−i (Ω)

i=1

F

=


IN 2 (:, Ω̄) A1 (:, Ω̄)

...


Ap (:, Ω̄)

(22)

12

are equal, which is added as constraint to the least squares
(22). The estimation of the left and right factor matrices along
with the missing data still belongs to the class of multi-convex
problems. The alternating minimization algorithm consists of
minimizing each set of variables while the others are fixed and
iterate until convergence to a stationary point (that might be
different from the global minimum). The least-squares (22)
is solved after line 19 of Algorithm 1. The computational
complexity of this step depends however on the amount
of missing data. The number of unknowns is of the order
Nt αN 2 with α the ratio of missing entries. The quality of
the prediction as a function of α is illustrated for a practical
scenario in Section VII. Theoretical bounds on α are not
derived in this paper.

unique as highlighted with Figure 3 and Figure 4, and both
solutions are equivalent as they provide a similar predictionerror (up to machine precision). The case (p, r) = (1, 2) is
analyzed in Figure 5 and Figure 6. We observe in the two latter
figures that using normalization affects the convergence speed
to the global minimum. In this example, about 500 iterations
were required: this observation is very much case dependent.
However, for all experiments carried out, the non-normalized
algorithm converges to the global minimum in few iterations.
10 0
No normalization
Correct normalization

Cost function residual

There is redundancy in the coefficients xk , i.e the block antidiagonals of the matrix


sp+1 (Ω̄) . . .
sNt (Ω̄)

 

..
..
..
xp . . . xNt = 

.
.
.
s1 (Ω̄)
. . . sNt −p−1 (Ω̄)

10 -5

10 -10

10 -15
0

2

4

6

8

10

12

14

16

18

20

Iteration

The proposed QUARKS identification method is first illustrated with a randomly generated VARX model and second,
with an Adaptive Optics (AO) application. The experiments
are carried out on MatlabR2016b in all subsections where a
scalability analysis is performed on a desktop computer with
a CPU Intel Xeon E5-2609V2/2.5 GHz with 8GB of RAM.
A. Case study 1: Randomly generated VARX model
We first illustrate the convergence of Algorithm 1 with
different normalizations for a randomly generated QUARKS
model whose temporal order and Kronecker rank is known.
The model structure is the following:
Sk =

p X
r
X

Figure 3. Evolution of the cost function as a function of the number
of iterations with 2 normalizations (no normalization, normalization as in
Algorithm 1). The pair (p, r) is set to (2, 1). The maximum number of ALS
iterations is set to 20.

10 0

ca k22
Residual of the fixed point: kMa − M

VII. N UMERICAL EXAMPLES

No normalization
Correct normalization

10 -5

10 -10

10 -15

(j)
(j)
M(ai )Uk−i M(bi )

0

2

4

6

8

10

12

14

16

18

20

Iteration

i=1 j=1
(j)

(j)

where N = 10, Sk ∈ R10×10 . The factor matrices M(ai )
(j)
and M(bi ) are generated with a random Toeplitz pattern and
such that its entries decay away from the diagonal. The input
is a white Gaussian noise. The number of temporal samples Nt
is set to 100 × N pr. Two scenarios were tested to analyze the
influence of the temporal and spatial order (p, r). In Figure 3
and Figure 4, the pair (p, r) is set to (2, 1). Figure 3 plots
the residual of the QUARKS cost function as a function of
the iteration number for both normalized and non-normalized
algorithms. Convergence to the global minimum is observed
for both cases. The convergence towards a unique fixed point
is shown with Figure 4 that displays the least-squares residual
d
(j)
(j)
between the true value M(ai ) and the estimate M(ai ) (as
defined in lines 29-30 of Algorithm 1). When normalizing
d
(j)
the columns of the factor matrix, M(ai ) converges to the
true value while it is not the case for the non-normalized
version. Although both algorithms reach the global minimum,
the solution to the QUARKS identification problem is not

Figure 4. Evolution of the least squares between the true value M(ai )
d
(j)
and the estimate M(ai ) (as defined in lines 29-30 of Algorithm 1) as a
function of the number of iterations with 2 normalizations (no normalization,
normalization as in Algorithm 1). The pair (p, r) is set to (2, 1). The
maximum number of ALS iterations is set to 20.

The algorithm converges in a monotonous manner to a fixed
point for both pairs (p, r). In the next example, we illustrate
the performance on a practical case study for different scenarii
of normalization and regularization.
B. Case study 2: Adaptive optics
When a light beam with an originally flat wavefront passes
through a turbulent medium, the wavefront gets distorted both
in time and space. AO is the system that corrects in real-time
atmospheric aberrations to obtain high resolution images. The
incoming distorted wavefront is flattened into a wavefront with
minimal aberrations using a deformable mirror operation in
an H2 -optimal feedback loop, [41]. However, the turbulence
evolves over time and accurate wavefront predictions are

13

Cost function residual

10 0

No normalization
Correct normalization

10 -5

10 -10

10 -15
0

100

200

300

400

500

600

700

Iteration

Figure 5. Evolution of the cost function as a function of the number
of iterations with 2 normalizations (no normalization, normalization as in
Algorithm 1). The pair (p, r) is set to (1, 2). The maximum number of ALS
iterations is set to 1000.

a grid of Nφ × Nφ samples and is denoted with φk when
lifted into a vector. In this paper, two layers of turbulence with
different statistics and windspeed are located on conjugated
planes and added up to form the wavefront measured by the
sensor. The atmospheric turbulence is a stochastic process,
therefore 50 realizations are carried out.
The wavefront cannot be measured itself, which is why a
wavefront sensor called Shack-Hartmann is considered here. It
consists of a N × N lenslet array parallel to a camera located
at the focal point of the lenslets. Each lenslet corresponds to
N −1
an identical number of wavefront points, nφ + 1 = φN . The
sensor measures locally the first derivative of the wavefront in
both horizontal and vertical direction. For one lenslet located
at the position (i, j), the noise-free measurement sk,(i,j) ∈ R2
is such that:
ZZ
I
sk,(i,j) =
∇φk dΣi,j =
φk dr
Σi,j

ca k22
Residual of the fixed point: kMa − M

10 0

10

-5

10 -10

No normalization
Correct normalization

10 -15
0

100

200

300

400

500

600

700

Iteration

Figure 6. Evolution of the cost function as a function of the number
of iterations with 2 normalizations (no normalization, normalization as in
Algorithm 1). The pair (p, r) is set to (1, 2). The maximum number of ALS
iterations is set to 1000.

required. A number of predictive methods have been designed
to compensate for turbulence, and each of them establish a
different trade-off between accuracy in minimizing the residual
phase variance and computational complexity. A state-space
algorithm has been derived in [41] to identify the turbulence
dynamics based on the wavefront sensor measurements. It
however does not scale to large-scale systems and alternatives
have been proposed using e.g diagonal VAR models, [39]. The
spatial covariance of the wavefront is not sparse, and hence
the latter modeling ignores the spatio-temporal correlation.
A data-sparse alternative for modeling large-scale turbulent
phase screens from sensor data using QUARKS identification
is proposed as a numerical illustration of a practical scenario.
The atmospheric turbulence is a stochastic and stable process that has been modeled via state-space [38] or using multiscale PCA for generating large-scale phase screens, [42]. Such
algorithms rely on the frozen-flow hypothesis, which states
that the time it takes for the turbulence layer to cross the line of
sight of the telescope is significantly larger than the sampling
frequency. The phase deformations are here generated accordingh to [42]
i using the
h knowni wavefront covariance matrices
T
E φk φk and E φk φTk+1 based on von Karman firstprinciples models. More specifically, only the low resolution
process is used based on the Fast Fourier Transform Moving
Average generator. The wavefront is spatially discretized on

∂Σi,j

where Σi,j is the surface over the lenslet (i, j). The wavefront
is defined on a discrete high-resolution grid, hence, the slope
measurements sk,(i,j) are computed with:


sx,k,(i,j)
sk,(i,j) :=
sy,k,(i,j)
Pnφ
Pnφ

xj =1 φk,(i,j),(nφ ,xj ) − Pxj =1 φk,(i,j),(1,xj )
P
=
nφ
nφ
xi =1 φk,(i,j),(xi ,nφ ) −
xi =1 φk,(i,j),(xi ,1)
where φk,(i,j) ∈ Rnφ ×nφ denotes the local wavefront measured by the lenslet (i, j). Photon noise and readout noise
are modeled with a zero-mean white Gaussian noise and
covariance matrix σv2 I2N 2 .
Unless mentioned otherwise in the next paragraphs, the
default parameters for AO simulations are listed in Table I.

Table I
D EFAULT PARAMETERS FOR AO SIMULATION .
Turbulence
Number of layers
Fried parameter, r0
Outer scale, L0
Number of phase points per lenslet, nφ
Horizontal windspeed,
Number of realizations,
Telescope
Telescope aperture, D
Sampling frequency
Number of lenslets, N
Number of sensor measurements, 2N 2
Signal-to-Noise Ratio, SNR

2
{0.2, 0.4} [m]
10 [m]
3
{1, 2} [points/sample]
50
1 [m]
500 [Hz]
10
200
15 [dB]

C. Benchmark methods and quality criteria
Three methods for identification are compared:
1) unstructured least squares
min
Ai

Nt
X
k=p+1

ksk −

p
X
i=1

Ai sk−i k22

(23)

14

2) regularized sparse least-squares using [40]:
min
Ai

Nt
X
k=p+1

ksk −

p
X

Ai sk−i k22 + τ

i=1

p
X

kvec(Ai )k1

i=1

(24)
where τ is a regularization parameter.
3) QUARKS identification (15) with Algorithm 1 without
normalization. No knowledge of the normalization coefficients is available. Algorithm 1 is initialized only
once, randomly. The stopping criterion parameters  and
Imax are set respectively to 10−5 and 3. The maximum
number of iterations `max is 100. In Subsection E,
the hyperparameters were randomly searched within the
bounds mentioned in Section III and within the range
[0, 5] for (λ, µ): the set of hyperparameters over 20
realizations that yields the lowest prediction-error is
selected. The curse of dimensionality that appears when
choosing hyperparameters with grid search is bypassed
with random search, [43]. Bayesian optimization or
online non-linear optimization for hyperparameter estimation are outside the scope of this paper.
The performances are checked on a validation dataset containing 5×103 temporal points. The results are discussed based
on the Variance Accounted
Pp c For (VAF) between the signals
sk+1 and b
sk+1 = i=1 A
i sk−i :
PNt
1

 
sk k22 
k=1 ksk − b
Nt
×
100
VAF(sk , b
sk ) = max 0, 1 −
P
Nt
1
2
k=1 ksk k2
Nt
The VAF between two identical signals sk and b
sk reaches
100%.
In Subsection D, we compare three structures for large-scale
VARX modeling and analyze the impact of increasing the spatial order r for QUARKS identification. In Subsection E, we
investigate the impact of temporal and spatial regularization
for varying SNR conditions. We illustrate the computational
complexity analysis from V.B with timing experiments for a
range of network sizes in Subsection F. Last, the missing data
algorithm as explained in VI.B is studied and its robustness
with respect to the ratio of unknown nodes is analyzed.

truncation on the elements of the factor matrices is done for the
Kronecker model. For example, a total of 500 non-zero values
are necessary to build the Kronecker factors associated to A1
and reaches 85.54% accuracy. The VAF obtained with the
sparse identification decreases with increasing regularization
parameter τ as expected while the number of non-zero entries
decreases for a high prior on sparsity. This trade-off between
the complexity of the model and the accuracy of the prediction
error is present in the QUARKS modeling as well. While
the estimated matrix with `1 minimization tries to reduce
the number of non-zero entries, the matrix obtained with
QUARKS modeling does not exhibit sparse patterns but a
prominent multi-level structure. The lower the spatial order r,
the lower the model complexity and the higher the prediction
error is.
90
r=2

88

r=3

r=1

r=4
r=5

86
84

=1

= 6.31
= 251.2

= 39.8

82
80
78
76
74
72
70
10 3

10 4

10 5

Figure 7. Variance Accounted For (%) versus complexity of model. A blue
cross corresponds to an estimate with given Kronecker rank. Each red cross
corresponds to a regularization parameter τ on the sparsity prior in (24). Two
points are not visible on the plot: (τ, %non-zero values, VAF) ∈ {(1.5849 ×
103 , 389, 45.3), (104 , 124, 0)}.

Figure 8 displays the coefficient-matrix A1 along with the
(1)
(1)
factor matrices M(a1 ) and M(b1 ), which gives insight on
the turbulence spatial correlations.
-1.5

D. Illustration of QUARKS identification

-1

200

-2

400

1

5
-2

10

5
0

15

-2.5

20

-3

-3

3

The identification set contains 5 × 10 temporal measurements. The temporal order of the VAR model is set to 2.
We first choose a Kronecker rank within r = {1, ..., 5}. The
parameters λ, µ in (15) are set to 0. The minimization (24)
is solved for τ = logspace(0, 4, 8). We define a measure that
we call model complexity as the number of non-zero entries
needed to construct the p coefficient matrices. For example,
the complexity of a QUARKS model is at most 2prN 2 (only
the non-zero elements of the factor matrices) while it reaches
a total of pN 4 for the full least squares estimation. It is
illustrated in Figure 7 that displays the VAF with respect
to the number of non-zero elements (with truncated entries
at 1% of the maximum value) needed to construct the full
coefficient matrix A1 . The prediction error is computed on
a validation dataset after truncation. We emphasize that no

600

-4

800

-5

10
-1

25
-3.5

15

30
-6

1000

-2

-4

35
-7

40

-4.5

-8

45

-5

1400

-9

50

1600

-10

55

20
-3

1200

25
-4

-5.5
-6

200

400

600

800

1000 1200 1400 1600

10

20

30

40

50

5

10

15

20

25

Figure 8. Identification of a QUARKS model with (p, r, N ) = (2, 2, 29)
and 1682 output channels. The absolute value of the entries in A1 (left), in
(1)
(1)
M(a1 ) (middle), M(b1 ) (right) are displayed in log10.

E. Influence of the hyperparameters
The regularization with rs and rt in (13) is the most
beneficial with short data batches or in noisy environments.
The difference with the case (λ, µ) = (0, 0) is all the
t
more significant when the ratio NNrp
is approximately 1. The
parameters for this subsection are gathered in Table II.

15

Table II
PARAMETERS FOR VII.E.
Sensor
Signal-to-Noise Ratio, SNR [dB]
QUARKS
Temporal order, p
Spatial order, r
Number of points for identification, Nt
Number of Monte-Carlo simulations

1.5
1

{5, 10, 15, 20, 25, 30, 35, 40}
0.5

4
2
500
50

0
-0.5
-1

100

-1.5

90

-2
80

-2.5
0.6

VAF(%)

70

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

60

Figure 10. Evolution of the computational time with respect to the size
of the 2D array. The linear model fitted with the QUARKS method is:
log10(Time) = 2.55 × log10(N ) − 3.10, σ = 0.34 while it is:
log10(Time) = 5.03 × log10(N ) − 5.68, σ = 0.27 with the unstructured
least-squares (LS).

50
40
30
20
5

10

15

20

25

30

35

40

45

50

SNR(dB)

Figure 9. Variance Accounted For (%) versus the signal-to-noise ratio. Red:
without regularization nor normalization. Blue: with both regularization and
normalization to 1.

Figure 9 displays the VAF on validation data with and
without regularization. Regularizing the cost function in noisy
situations and with relatively few data samples leads to
substantial improvements over the non-regularized QUARKS
identification. It especially reduces the variance of the prediction error while the performance of the non-regularized version
with few temporal samples is very unreliable. Random search
has interesting performances as it exploits the fact that some
hyperparameters may not contribute a lot for obtaining good
solutions in the example at hand.

G. Missing data
Last, we analyze the missing data algorithm as described
in Subsection VI.B. We choose (p, r) = (2, 2). Although
the SNR in this paragraph is the default one, i.e 15dB,
the number of identification points is 5 × 103 and hence
much larger than in Subsection E which allows for a much
better identification. The regularization parameter β range in
{0, logspace(−2, 2, 5)}. For a given experiment, the position
of the missing entries are fixed for all temporal points in the
identification batch. 70 experiments are carried out and each
of them has a different pattern for the position of missing data.
The VAF is plotted as a function of the percentage of missing
entries in Figure 11. The prediction performances of the
100
90
80
70
60

VAF(%)

F. Scalability
One advantage of the new modeling paradigm is to reduce
the computational complexity for estimating large-scale VARX
models. No regularization is considered in this section in order
to analyze whether the QUARKS identification in Algorithm
1 scales with O(N 3 Nt ) as highlighted in V.B. The temporal
order is set to 4 and Kronecker rank to 2. The number of
lenslets N belongs to the range [5 : 2 : 29], which implies 2 ×
{52 , ..., 292 } sensor signals at each time sample. The number
of time samples for QUARKS identification is such that Nt =
10prN while it is Nt = 50N 2 . These values were fixed such
that the prediction-error is similar for both methods. The linear
model fitted in Figure 10 for the QUARKS has a regression
coefficient of 3.27 (with standard deviation σ = 0.51) while
the unstructured estimation has a coefficient of 5.18 (σ =
0.50). Although the QUARKS implementation includes many
loops and would take advantage of a C implementation, the
reduction in the regression coefficient is significant using the
Kronecker-based identification.

50
40
30
20
10
0
0

10

20

30

40

50

60

70

80

90

100

Percentage of missing entries in Sk

Figure 11. Variance Accounted For for the following methods: (Red) nonregularized VAR-2 with unstructured coefficient-matrices (VAR2-p = 2),
(Green) QUARKS with β = 0 and no normalization, (Blue) QUARKS with
optimal β and no normalization. For each experiment, there is a pattern of
missing entries which gives rise to one pair of blue and green crosses.

alternating minimization is poor when β = 0 as it decreases
very sharply to 0 when the ratio of missing entries increases.

16

However, for an optimal β, the VAF raised significantly and
reaches better prediction errors than with the unstructured
unregularized VAR estimation. Nonetheless, it should be kept
in mind the cost of computing the least-squares (22) when the
ratio of missing data increases hampers the scalability of the
proposed method. A remedy to such a situation was studied
in [24] in the context of blind identification.
VIII. C ONCLUSION
In this paper, the class of Kronecker networks is defined
and the VAR modeling part is investigated. Each coefficientmatrix of the VAR model is approximated with a sum of
few Kronecker matrices which offers high data compression
for large networks. Estimating in least-squares sense the data
matrices give rise to a bilinear problem which is addressed
using Alternating Least Squares (ALS). The convergence of
ALS to a global minimum was proven assuming persistency of
excitation and rank conditions on the initial guesses. The ALS
framework provides an elegant way of dealing with missing
sensor data and adding further information on the factor
matrices, e.g Toeplitz, banded, sparse. Numerical examples
on atmospheric turbulence prediction, e.g for AO control,
demonstrates the high compression capabilities of this model
as well as its scalability for larger networks.
The algorithm has been presented for 2D dynamical systems
and can easily be generalized to higher dimensions by using
a Kronecker product of multiple matrices instead of only two
matrices in which case larger compression rates are achieved.
Such higher order modeling for 2D arrays is obtained by
tensorizing the sensor data sk and allows to establish a new
trade-off between accuracy and computational complexity.
The new paradigm of Kronecker networks identification
opens up wide research possibilities in modeling large-scale
dynamical systems. Future research consists in expanding
on the QUARKS to identify V-ARX, V-ARMAX and BoxJenkins as well as subspace identification models. The latter
is investigated in the companion paper [25].

In this Appendix, we derive the proof of convergence for
the regularized Alternating Least Squares in matrix form. The
proof builds on [36] and therefore we only highlight the main
changes here compared to the vector form. First, an inner
3
product for matrices in R(N (Nt −1)+N )×N is defined.
Definition A1. Let X, Y ∈ RN ×N and denote their columns
with xi , yi . For two matrices X̄, Ȳ such that


 IN ⊗ x1

I ⊗ Ũ


..
X̄ = √ N


.
µ(IN ⊗ IN 2 )
IN ⊗ x N

hX̄, Ȳi = λmax (ŨT Ũ + µIN 2 )vec(X)T vec(Y)

3

Lemma A1. For the matrix X̄ ∈ R(N (Nt −1)+N )×Npand the
inner product in Definition 4, the quantity kX̄k2 = hX̄, X̄i
3
is a norm on R(N (Nt −1)+N )×N .
Proof. The proof contains 4 points.
1) Positiveness: kX̄k2 is positive because λmax and the
2-norm for vectors k.k2 are both positive.
2) If kX̄k2 = 0, and with ŨT Ũ+µIN 2 6= 0, then kxk2 = 0
and x = 0. This implies that X̄ = 0.
3) Let α ∈ R.


kαX̄k22 = λmax α2 (ŨT Ũ + µIN 2 ) kxk22
= |α|kXk22
4) Triangular inequality. We can write:
q
kX̄ + Ȳk2 =
λmax (ŨT Ũ + µIN 2 )kx + yk2
q
λmax (ŨT Ũ + µIN 2 )(kxk2 + kyk2 )
≤
≤

kX̄k2 + kȲk2

k.k2 is therefore a norm.
.
Recall the definition of the two sets, associated to each of
the variables b and a:
Xb

2
= {b
b ∈ RN |∀i ∈ {1, ..., N }, kbbi k2 = kbi k2 , bb1 > 0}

Xa

= {b
a ∈ RN |∀i ∈ {1, ..., N }, kabi k2 ≤ kai k2 }

2

Let a ∈ Xa , b ∈ Xb . The noise term Ṽ is defined similarly
as S̃ from the noise components ek . We denote ỹ := vec(Ỹ).
A. F is an operator mapping Xa to Xa
We write the solution of one least-squares update with:
b
F3 (b)


¯ T M̄
¯ )−1 M̄
¯ T ỹ
=
IN ⊗ (M̄
b̂
b̂
b̂
 


T ¯
Ṽ
−1 ¯ T
¯
¯
=
IN ⊗ (M̄b̂ M̄b̂ ) M̄b̂
(IN ⊗ M̄b )a + vec
0

b =
a

A PPENDIX A
P ROOF FOR T HEOREM 1

and similarly for Ȳ, the inner product on R(N (Nt −1)+N
is defined with:

¯ has the structure of X̄ in the
For example, the matrix M̄
a
above definition. Denote:


I ⊗ Ũ
M= √ N
µ(IN ⊗ IN 2 )

3

)×N

or equivalently:
¯ a +v
¯ T M̄
¯ )−1 M̄
¯ T M̄
b i = (M̄
a
b i
i
b̂
b̂
b̂
It comes,
lim kb
ai k2

Nt →∞

≤

¯ T M̄
¯ )−1 M̄
¯ T M̄
¯ k ka k
k(M̄
b 2
i 2
b̂
b̂
b̂
Nt →∞

T
T
¯ M̄
¯ )−1 M̄
¯ k kv k
+k(M̄
2
i 2
b̂
b̂
b̂
lim



Moreover:
¯ T M̄
¯ )−1 M̄
¯ Tk
lim k(M̄
2
b̂
b̂
b̂

Nt →∞

≤

¯ k
kM̄
b̂ 2
Nt →∞ ¯ T ¯
kM̄ M̄ k
lim

b̂

b̂ 2

17

C. Upper bound on Q

Using the definition of the inner product,

b (`+1) =
We now introduce the quantity Q = k dF
db
a k2 . From a
F3 (F2 (F1 (b
a(`) ))), we decompose:

¯ k
kM̄
1
b̂ 2
q
=
T
¯ M̄
¯ k
kM̄
kb
bk2 λmax (ŨT Ũ + µIN 2 )
b̂
b̂ 2
We assumed the noise has a finite variance when Nt goes to
infinity, which implies:

≤

T

T

Q =

¯ k kv k = 0
¯ M̄
¯ )−1 M̄
lim k(M̄
2
i 2
b̂
b̂
b̂

Nt →∞

Therefore,
lim kb
ai k2

Nt →∞

≤

¯ T M̄
¯ k ka k
¯ T M̄
¯ )−1 M̄
lim k(M̄
b 2
i 2
b̂
b̂
b̂

Nt →∞

¯ T M̄
¯ k
kM̄
b 2
b̂
≤
lim
¯ k kai k2
Nt →∞ kM̄ M̄
b̂

≤

b̂ 2

b T bk2
kb
ka k
lim
b T bk
b 2 i 2
Nt →∞ kb

The inequality is strict if b
b 6= b.
We have b ∈ Xb , which implies kb
bk2 = kbk2 and therefore,
using Lemma A2, kb
ai k2 ≤ kai k2 when Nt goes to infinity.
Then, the functional F(.) maps Xa to Xa .
d
B. Proof that kb
ai k2 = kai k2
ni k2 = kbi k2 and kb
For Nt going to ∞:


¯ T M̄
¯ )−1 M̄
¯ T M̄
¯ a
b = IN ⊗ (M̄
a
b̂
b̂
b̂
b̂

(25)

¯ on both left sides and using similar
Multiplying by IN ⊗ M̄
b̂
¯ M̄ . The right-hand side
¯ M̄ = M̄
arguments as in [36], M̄
â
b
a
b̂
term reads:


IN ⊗ b1


..
¯ M̄ = M 
M̄

 a1 . . . aN
b
a
.
IN ⊗ bN
and hence, for all i, j:
M(IN ⊗ bi )aj = M(IN ⊗ bbi )abj
The matrix M is full column rank, it follows:
(IN ⊗ bi )aj
bi aj k

(IN ⊗ bbi )abj
= bbi abj k

=

Therefore, since aj k ∈ R and b ∈ Xb , it follows:
kbi k2 |aj k | = kbbi k2 |abj k | and then, |aj k | = |abj k |, for all k in
the set {1, ..., N }. Finally, it comes kaj k2 = kabj k2 . A similar
reasoning starting from the relation between b
b and b yields
cni k2 = kbi k2 .
kb

2

2

dF1
db
a

2

Lemma A3. Let f (.) be defined with: f (b
b) := IN ⊗
¯ T M̄
¯ )−1 M̄
¯ T . Under Assumption A2, the magnitude of
(M̄
b̂
b̂
b̂
the directional derivative of f (b
b) along a vector u attains its
maximum when u is in the same direction as b
b.
Hence, similarly as in [36],
→
−
−
k→
a T ak2
k b T bk2
kak
kbk2
Q ≤
2
T
b k2
kb
aT a
kb
b b
bk2
≤

We conclude with the lemma:
b ∈ RN . If kb
Lemma A2. Let  ∈ {−1, 1} and b, b
bk2 = kbk2
b
and b 6= b, then
T
T
kb
b bk2 ≤ kb
b b
bk2

dF db
b db
bn
b
b
db
db dbn a
dF3
dF2
b
db 2 db
bn

T
aT ak2
kb
b bk2 kb
T
b k2
aT a
kb
b b
bk2 kb

(26)

obtained using b ∈ Xb , a ∈ Xa . Using Lemma A2, we
conclude that Q < 1. Hence, F is a contraction map on Xa
and the ALS algorithm converges to the fixed point.
R EFERENCES
[1] R. Gilmozzi, J. Spyromilio, “The European Extremely Large Telescope
(E-ELT)”, The Messenger 12 7, 11, 2007.
[2] G. Derado, F. D. Bowmqn, C.D. Kilts,“Modeling the spatial and temporal
dependence in fMRI data”, Biometrics, vol. 66, no. 3, pp. 949-957, 2010.
[3] C. Kulcsár, H-F Raynaud, C. Petit, J-M Conan, “Minimum variance
prediction and control for adaptive optics”, Automatica 48 (2012)
1939–1954.
[4] R. P. Roesser, “A Discrete State-Space Model for Linear Image Processing”, IEEE Trans. on Automatic Control, Vol. AC-20, No. 1, 1975.
[5] E. Fornasini, G. Marchesini, “Doubly-Indexed Dynamical Systems: StateSpace Models and Structural Properties”, Math. Systems Theory 12, 5972, 1978.
[6] A. Alenany, G. Mercère, J. A. Ramos, “Subspace Identification of 2D CRSD Roesser Models With Deterministic-Stochastic Inputs: A State
Computation Approach”, IEEE Trans. on Control Systems Technology,
2016.
[7] A. Chiuso, G. Pillonetto, “A Bayesian approach to sparse dynamic
network identification”, Automatica 48 (2012) 1553-1565.
[8] J. K. Rice, “Efficient Algorithms for Distributed Control: A Structured
Matrix Approach”, PhD thesis, Delft Center for Systems and Control, TU
Delft, 2010.
[9] C. Yu, M. Verhaegen. “Local subspace identification of distributed homogeneous systems with general interconnection patterns”, IFAC symposium
SYSID, Beijing (pp. 585–589), 2015.
[10] M. Ali, S. S. Chughtai, H. Werner, “Identification of spatially interconnected systems”, in Proceedings of the 48ht IEEE Conference on Decision
and Control, pp. 7163-7168., 2009.
[11] P. Massioni, M. Verhaegen, “Subspace identification of circulant systems”, Automatica 44 (2008) 2825–2833.
[12] P. Massioni, “Distributed control for alpha-heterogenous dynamically
coupled systems”, Systems & Control Letters, vol. 72, pp 30–35, 2014.
[13] C. F. van Loan, N.P. Pitsianis, “Approximation with Kronecker products”, Linear Algebra for Large Scale and Real Time Applications,
Kluwer Publications, Dordrecht, vol.25, pp 293–314, 1992.
[14] C. F. van Loan, “The ubiquitous Kronecker product”, Journal of Computational and Applied Mathematics, vol.123, pp 85–100, 2000.
[15] L. Grasedyck, D. Kressner, C. Tobler, “A literature survey of low-rank
tensor approximation techniques”, GAMM-Mitteilungen, vol. 36, no. 1,
pp. 53-78, 2013.
[16] F. Bijma, J. de Munck, R. Heethaar, “The spatiotemporal MEG covariance matrix modeled as a sum of Kronecker products”, NeuroImage, vol.
27, pp. 402-415, 2005.

18

[17] T. Tsiligkaridis, A.O. Hero, “Covariance Estimation in High Dimensions
Via Kronecker Product Expansions”, IEEE Trans. on Signal Processing,
Vol. 61, No. 21, 2013.
[18] C. Hansen, J. G. Nagy and D. P. O’Leary, Deblurring images : matrices,
spectra, and filtering. Fundamentals of Algorithms 3, SIAM, Philadelphia,
2006.
[19] A. Doostan, G. Iaccarino, “A least-squares approximation of partial
differential equations with high-dimensional random inputs”, Journal of
Computational Physics 228(2009) 4332-4345.
[20] L. Giraldi, A. Nouy, and G. Legrain “Low-rank approximate inverse
for preconditioning tensor-structured linear systems”, SIAM Journal on
Scientific Computing, vol. 36, no. 4, pp. 1850-1870, 2014.
[21] G. Beylkin and M. J. Mohlenkamp, “Algorithms for numerical analysis
in high dimensions”, SIAM Journal on Scientific Computing, vol. 26, no.
6, pp. 2133-2159, 2005.
[22] A. Cichocki et al., “Tensor Networks for Dimensionality Reduction and
Large-scale Optimization: Part 2 Applications and Future Perspectives,
Foundations and Trends in Machine Learning, Vol 9, Issue 6, May 2017.
[23] P. D. Hoff, “Multilinear tensor regression for longitudinal relational
data”, The Annals of Applied Statistics, Vol. 9, No. 3, 1169-1193, 2015.
[24] M. Boussé, O. Debals, L. De Lathauwer, “A Tensor-Based Method
for Large-Scale Blind Source Separation Using Segmentation”, IEEE
Transactions on Signal Processing, Vol. 65, No. 2, 2017.
[25] B. Sinquin, M. Verhaegen, “K4SID: Large-Scale Subspace Identification
with Kronecker modeling”, Submitted to IEEE Trans. on Automatic
Control, 2017.
[26] B. Sinquin, M. Verhaegen, “Kronecker-ARX Models in Identifying (2D)
Spatial-Temporal Systems”, IFAC World Congress, Toulouse, France,
2017.
[27] T. Hastie, R. Tibshirani, J. Friedman, “The Elements of Statistical
Learning”, Springer, 2008.
[28] T. Chen, H. Ohlsson, L. Ljung, “On the estimation of transfer functions,
regularizations anf Gaussian processes-Revisited”, Automatica 48 (2012)
1525-1535.
[29] G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, L. Ljung, “Kernel
methods in system identification, machine learning and function estimation: A survey”, Automatica 50 (2014) 657-682. Automatica 72 (2016)
217-229.
[30] M. Udell, C. Horn, R. Zadeh, S. Boyd, “Generalized Low-Rank models”,
Foundations and Trends in Machine Learning, Vol. 9, No. 1 (2016), 1118.
[31] P. Baldi, K. Hornik, “Neural networks and Principal Components Analysis: Learning from Examples Without Local Minima”, Neural Networks,
Vol 2, pp. 53-58, 1989.
[32] C. E. Rasmussen, C. K. I. Williams, “Gaussian Processes for Machine
Learning”, the MIT Press, ISBN 026218253X, 2006.
[33] F. Ding, T. W. Chen, “Iterative least-squares solutions of coupled matrix
equations”, Systems & Control Letters, 54 (2005) 95-107.
[34] J. Bruls, C. T. Chou, B. R. J. Haverkamp, M. Verhaegen, “Linear and
nonlinear system identification using separable least-squares”, European
Journal of Control, 5(1):116-128, 1999.
[35] N. Li, S. Kindermann, C. Navasca, “Some convergence results on the
Regularized Alternating Least-Squares method for tensor decomposition”,
Linear Algebra and its Applications 438 (2013) 796-812.
[36] G. Li, G. Wen, A. Zhang, “Fixed point iteration in identifying bilinear
models”, Systems& Control Letters 83(2015) 28-37.
[37] A. Granas, J. Dugundji Fixed Point Theory. Springer-Verlag, New York,
2001.
[38] A. Beghi, A. Cenedese, A. Masiero, Stochastic realization approach to
the efficient simulation of phase screens. J. Opt. Soc. Am. A, volume 25,
pages 515–525, 2008.
[39] C. Correia, K. Jackson, J-P. Veran, D. Andersen, O. Lardière, C. Bradley
Static and predictive tomographic reconstruction for wide-field multiobject adaptive optics systems. J. Opt. Soc. Am. A, vol.31, 2014.
[40] S.-J. Kim, K. Koh, M. Lustig, S. Boyd, D. Gorinevsky An Interior-Point
Method for Large-Scale l1-Regularized Least Squares. IEEE Journal on
Selected Topics in Signal Processing, 1(4):606-617, December 2007.
[41] K. Hinnen, M. Verhaegen, N. Doelman, “A Data-Driven H2 -Optimal
Control Approach for Adaptive Optics”, IEEE Trans. on Control System
Technology, Vol. 16, No. 3, 2008.
[42] A. Beghi, A. Cenedese, and A. Masiero, “Multiscale stochastic approach for phase screens synthesis,” Applied Optics, vol. 50, no. 21, pp.
4124–4133, 2011.
[43] J. Bergstra, Y. Bengio, “Random Search for Hyper-Parameter Optimization”, Journal of Machine Learning Research 13(2012) 281-305.

Baptiste Sinquin received his engineering degree from Ecole Centrale Lyon,
France, in 2014. He is currently pursuing the PhD degree at the Delft
Center for Systems and Control, TU Delft, the Netherlands, where he studies
system identification of large-scale systems. His research activities include
the validation of the identification and control methods in a laboratory
environment dedicated for large-scale adaptive optics.

Michel Verhaegen received an engineering degree in aeronautics from the
Delft University of Technology, Delft, The Netherlands, in 1982 and the
doctoral degree in applied sciences from the Catholic University Leuven,
Belgium, in 1985. From 1985 to 1994, he was a Research Fellow of the
U.S. National Research Council (NRC), affiliated with the NASA Ames
Research Center in California, and of the Dutch Academy of Arts and
Sciences, affiliated with the Network Theory Group of the Delft University
of Technology. From 1994 to 1999, he was an Associate Professor of the
Control Laboratory, Delft University of Technology and was appointed as
Full Professor with the Faculty of Applied Physics, University of Twente,
The Netherlands, in 1999.
In 2001, he moved back to the University of Delft and is now a member
of the Delft Center for Systems and Control. He has held short sabbatical
leaves at the University of Uppsala, McGill, Lund, and the German Aerospace
Research Center (DLR) in Munich and is participating in several European
Research Networks. His main research interest is the interdisciplinary domain
of numerical algebra and system theory. In this field he has published over
130 papers. Current activities focus on the transfer of knowledge about
new identification, fault tolerant control and data-driven controller design
methodologies to industry. Application areas include adaptive optics.

