Tensor Network alternating linear scheme for MIMO Volterra system
identification
Kim Batselier a , Zhongming Chen a , Ngai Wong a ,

arXiv:1607.00127v2 [cs.NA] 18 Oct 2016

a

The Department of Electrical and Electronic Engineering, The University of Hong Kong

Abstract
This article introduces two Tensor Network-based iterative algorithms for the identification of high-order discrete-time nonlinear multipleinput multiple-output (MIMO) Volterra systems. The system identification problem is rewritten in terms of a Volterra tensor, which is
never explicitly constructed, thus avoiding the curse of dimensionality. It is shown how each iteration of the two identification algorithms
involves solving a linear system of low computational complexity. The proposed algorithms are guaranteed to monotonically converge and
numerical stability is ensured through the use of orthogonal matrix factorizations. The performance and accuracy of the two identification
algorithms are illustrated by numerical experiments, where accurate degree-10 MIMO Volterra models are identified in about 1 second in
Matlab on a standard desktop pc.
Key words: Volterra series; tensors; MIMO; identification methods; system identification

1

Introduction

Volterra series [26] have been extensively studied and applied in applications like speech modeling [17], loudspeaker
linearization [13], nonlinear control [8], active noise control [22], modeling of biological and physiological systems [16], nonlinear communication channel identification
and equalization [4], distortion analysis [24] and many others. Their applicability has been limited however to “weakly
nonlinear systems”, where the nonlinear effects play a nonnegligible role but are dominated by the linear terms. Such
limitation is not inherent to the Volterra series themselves,
as they can represent a wide range of nonlinear dynamical
systems, but is due to the exponentially growing number of
Volterra kernel coefficients as the degree increases. Indeed,
assuming a finite memory M , the dth-order response of
a discrete-time single-input single-output (SISO) Volterra
system is given by

yd (t) =

M
−1
X

hd (k1 , . . . , kd )

k1 ,...,kd =0

d
Y

u(t − ki ),

i=1

where y(t), u(t) are the scalar output and input at time t respectively and the dth-order Volterra kernel hd (k1 , . . . , kd )
Email addresses: kim.batselier@eee.hku.hk (Kim
Batselier), zmchen@eee.hku.hk (Zhongming Chen),
nwong@eee.hku.hk (Ngai Wong).

Preprint submitted to Automatica

is described by M d numbers. For a multiple-input multipleoutput (MIMO) Volterra system with p inputs the situation
gets even worse, where the dth-order Volterra kernel for one
particular output is characterized by (pM )d numbers. This
problem is also commonly known as the curse of dimensionality.
The paradigm used in this article to break this curse is to
trade storage for computation. This means that all the
Volterra coefficients are replaced by only a few numbers,
from which all Volterra coefficients can be computed. This
idea is not new, e.g. Volterra kernels have been expanded on
orthonormal basis functions in order to reduce their complexity [2,6]. Tensors (namely, multi-dimensional arrays that
are generalizations of matrices to higher orders) are also
suitable candidates for this purpose. In [7] both the canonical polyadic [3, 10] and Tucker tensor decompositions [23]
were used. The canonical polyadic decomposition can suffer from instability however, and the determination of its
rank is known to be a NP-hard problem [12], which can be
ill-posed [5]. The main disadvantage of the Tucker decomposition of the Volterra kernels, which is in fact an expansion onto a set of orthonormal basis functions, is that it still
suffers from an exponential complexity.
This motivates us to develop and introduce a new description
of discrete-time MIMO Volterra systems in terms of particular Tensor Networks (TN) [18]. For the particular case
of multiple-input-single-output (MISO) systems, these Tensor Networks will turn out to be Trains (TTs) [19]. A TN

19 October 2016

representation does not suffer from any instability or exponential complexity and can represent all Volterra kernels
combined by O(d(pM + 1)r2 ) elements, where r is a to-bedetermined number called the TN-rank. TNs were originally
developed in the physics community. Of particular importance is the Density Matrix Renormalization Group (DMRG)
algorithm [25], which is an iterative algorithm originally developed for the determination of the ground state of a entangled multi-body quantum system. Its applicability however
is not limited to problems in quantum physics, as demonstrated by recent interest in the scientific computing community [11, 20, 21].

i3
i2
i1

21
22
23
24

A d-way or dth-order tensor is denoted A ∈ Rn1 ×n2 ×···×nd
and hence each of its entries ai1 i2 ···id is determined by d
indices. The numbers n1 , n2 , . . . , nd are called the dimensions of the tensor. A tensor is cubical if all its dimensions
are equal. A cubical tensor is symmetric when all its entries
satisfy ai1 i2 ···id = aπ(i1 ,i2 ,...,id ) , where π(i1 , i2 , . . . , id ) is
any permutation of the indices. An example 3-way tensor
with dimensions 4, 3, 2 is shown in Fig. 1. For practical purposes, only real tensors are considered. We use boldface
capital calligraphic letters A, B, . . . to denote tensors, boldface capital letters A, B, . . . to denote matrices, boldface
letters a, b, . . . to denote vectors, and Roman letters a, b, . . .
to denote scalars. The transpose of a matrix A or vector a
are denoted by AT and aT , respectively. The unit matrix of
order n is denoted In .

(1) we derive a new description of discrete-time MIMO
Volterra systems using the TN format,
(2) we derive two iterative MIMO Volterra identification
algorithms that estimate all Volterra kernels in the TN
format from given input-output data.
In each step of the iterative identification a small linear system needs to be solved. The main computational tools are the
singular value decomposition (SVD) and the QR decomposition [9]. These orthogonal matrix factorizations ensure the
numerical stability of the methods [11]. The first identification method, which is called the Alternating Linear Scheme
(ALS) method, has the lowest computational complexity
but assumes that the TN-ranks rk ’s are fixed. The second
identification method, called the Modified Alternating Linear Scheme (MALS), removes this limitation and allows for
the adaptive updating of the TN-ranks during the iterations.
The TN-ranks are determined numerically by means of a
SVD. This is reminiscent of the determination of the order
of linear systems in subspace identification algorithms [14].
Monotonic convergence of both the ALS and MALS methods under certain conditions is discussed in [21].

We now give a brief description of some required tensor
operations and properties. The generalization of the matrixmatrix multiplication to tensors involves a multiplication of
a matrix with a d-way tensor along one of its d possible
modes (equiv. indices/axes).
Definition 2.1 (Tensor k-mode product)( [15, p. 460])
The k-mode product B = A ×k U of a tensor A ∈
Rn1 ×···×nk ×···×nd with a matrix U ∈ Rpk ×nk is defined by

The outline of this article is as follows. In Section 2 we
give a brief overview of important tensor concepts, operations and properties. The MIMO Volterra TN framework is
introduced in Section 3. The two iterative identification algorithms are derived in Section 4 and applied on two examples in Section 5. To our knowledge, this is the only time
where MIMO Volterra systems of degree 10 were identified
in about 1 second on a standard computer. Matlab/Octave
implementations of our algorithms are freely available from
https://github.com/kbatseli/MVMALS.

2.1

5 9
6 10
7 11
8 12

17
18
19
20

Fig. 1. A 4 × 3 × 2 tensor where each entry is determined by three
indices i1 , i2 , i3 .

In this article, we adopt the DMRG method in the TN format for the identification of MIMO discrete-time Volterra
systems. The contributions of this article are twofold:

2

1
2
3
4

13
14
15
16

bi1 ···ik−1 jik+1 ···id =

nk
X

ujik ai1 ···ik−1 ik ik+1 ···id ,

(1)

ik =1

and B ∈ Rn1 ×···×nk−1 ×pk ×nk+1 ×···×nd .
The following illustrative example rewrites the familiar matrix multiplication as a 1-mode and 2-mode product.
Example 1 (Matrix multiplication as mode products) For
matrices A, B, C with matching dimensions we have that

Preliminaries

A ×1 B ×2 C := B A C T .

Tensor basics

An interesting observation is that the definition of the kmode product also includes the multiplication of a tensor

Tensors in this article are multi-dimensional arrays that generalize the notions of vectors and matrices to higher orders.

2

A with d vectors. The following example highlights a very
important case.

n3

n1

Example 2 Consider a symmetric d-way tensor A with dimensions n and a vector x ∈ Rn . The multidimensional


n1

r2
r2

A

r1

(2)

A (3)
n3

n2

is the scalar

A xd := A ×1 xT ×2 xT ×3 · · · ×d xT ,

Fig. 2. The TT-cores of a 3-way tensor A are two matrices
A(1) , A(3) and a 3-way tensor A(2) .

(2)

Definition 2.4 (Vectorization)( [15, p. 460]) The vectorization of a tensor A, denoted vec(A), rearranges all its entries into one column vector.

which is obtained as a homogeneous polynomial of degree
d in the variables x1 , . . . , xn .

Example 4 For the tensor in Fig. 1, we have

The Kronecker product plays a crucial role in our description
of MIMO Volterra systems.

vec(A) = reshape(A, [24, 1]) =

Definition 2.2 (Kronecker product)( [15, p. 461]) If B ∈
Rm1 ×m2 and C ∈ Rn1 ×n2 , then their Kronecker product
B ⊗ C is an m1 × m2 block matrix whose (i3 , i4 )th block
is the n1 × n2 matrix bi3 i4 C




b11 · · · b1n1
b11 C · · · b1n1 C
 . .

 .
.
..
..
 .
..
..  ⊗ C =  ..
.
.
 .


bm1 1 · · · bm1 n1

A

=

A (1)

n2

T

contraction of A with x = x1 x2 · · · xn

r1

1 2 · · · 24

T

.

The importance of the vectorization lies in the following
equation


vec(A ×1 U1 ×2 · · · ×d Ud ) = (Ud ⊗ · · · ⊗ U1 )vec(A).
(5)


.


bm1 1 C · · · bm1 n1 C

Observe how the order in the Kronecker product is reversed
with respect to the ordering of the mode products. Equation
(5) allows us to rewrite (2) as

(3)
We use the notation x d := x ⊗ x ⊗ · · · ⊗ x for the d-times
repeated Kronecker product. The mixed-product property of
the Kronecker product states that if A, B, C and D are
matrices of such sizes that one can form the matrix products
AC and BD, then
(A ⊗ B) (C ⊗ D) = AC ⊗ BD.



A xd = vec(A)T x d ,

(6)

which tells us how the k-mode products of a tensor A with
a vector x can be computed. When A in (5) is also a matrix,
we then obtain the following useful property

(4)
vec(U1 A U2T ) = (U2 ⊗ U1 ) vec(A).

Definition 2.3 (Reshaping)( [15, p. 460]) Reshaping is another often used tensor operation. The most common reshaping is the matricization, which reorders the entries of
A into a matrix. We adopt the Matlab/Octave reshape operator “reshape(A, [n1 , n2 , n3 · · · ])”, which reshapes the
tensor A into a tensor with dimensions n1 , n2 , n3 , . . .. The
total number of elements of A must be the same as n1 ×
n2 × n3 · · · .

2.2

(7)

Tensor Train decomposition

The Tensor Network representation used in this article lies
conceptually very close to the Tensor Train decomposition [19]. We therefore first discuss the Tensor Train representation. A TT-decomposition represents a d-way tensor A
in terms of d 3-way tensors A(1) , . . . , A(d) , which are called
the TT-cores. The kth TT-core has dimensions rk−1 , nk , rk ,
where rk−1 , rk are called the TT-ranks. Specifically, each
entry of A ∈ Rn1 ×···×nd is determined by

Example 3 We illustrate the reshaping operator on the 4 ×
3 × 2 tensor of Fig. 1



1 5 9 13 17 21


2 6 10 14 18 22


reshape(A, [4, 6]) = 
.
3 7 11 15 19 23


4 8 12 16 20 24

(1)

(2)

(d)

ai1 i2 ···id = Ai1 Ai2 · · · Aid ,
(k)

(8)

where Aik is the rk−1 ×rk matrix obtained from specifying
ik . Since ai1 i2 ···id is a scalar, it immediately follows that
r0 = rd = 1. The TT-decomposition is illustrated for a 3way tensor in Fig. 2. Note that when all TT-ranks are equal to
r, then the storage of a cubical d-way tensor with dimensions
n in the TT format needs O(dnr2 ) elements. Small TT-ranks

Probably the most important reshaping of a tensor is the
vectorization.

3

therefore result in a significant reduction of required storage
cost. In order to describe MIMO Volterra systems we will
need to remove the constraint that r0 = 1 and therefore
obtain a slightly more general TN. Observe that the removal
of this constraint implies that the left-hand side of (8) will
then be a r0 × 1 vector. The following multidimensional
contraction in the TN format turns out to be very important
for Volterra systems.

Definition 3.1 (SISO Volterra tensor) Given a discrete-time
SISO Volterra system of degree d and memory M as described in (10), then the d-way cubical Volterra tensor V of
dimension M + 1 is defined by

Lemma 2.1 ( [19, p. 2309]) Given a cubical tensor A ∈
Rn×···×n in the TT format and a vector x ∈ Rn . The multidimensional contraction A xd can then be computed as


T
ut := 1 u(t) u(t − 1) · · · u(t − M + 1)
∈ RM +1 .

y(t) = V udt ,
where

Axd = (A(1) ×2 xT ) (A(2) ×2 xT ) · · · (A(d) ×2 xT ).
(9)

The extension of this definition to the MIMO case is
straightforward. Indeed, the ith output yi (t) is then a multivariate polynomial of degree d in u1 (t), . . . , u1 (t − M +
1), u2 (t), . . . , u2 (t − M + 1), . . . , up (t − M + 1). This implies that the vector ut in Definition 3.1 needs to be extended
with these additional inputs, which results in an increase
in the dimension of the corresponding Volterra tensor. In
addition, the order of the Volterra tensor is incremented to
accommodate for the multiple number of outputs.

Again, in the MIMO case this contraction will result in a
r0 × 1 vector where now each entry is the evaluation of
a homogeneous polynomial. The numerical stability of the
two identification algorithms described in this article relies
on the TN-cores being either left or right orthogonal.

Definition 3.2 (MIMO Volterra tensor) Given a discrete
time p-input l-output Volterra system of degree d and memory M , then the d + 1-way Volterra tensor V of dimensions
l × pM + 1 × · · · × pM + 1 is defined by

Definition 2.5 (Left orthogonal and right orthogonal TNcores)( [11, p. A689]) A TN-core A(i) is left orthogonal if
it can be reshaped into an ri−1 ni × ri matrix A for which
AT A = Iri



y1 (t)


y2 (t)


y(t) :=  .  = V ×2 uTt ×3 uTt · · · ×d+1 uTt , (11)
 .. 


yl (t)

applies. Similarly, a TN-core A(i) is right orthogonal if it
can be reshaped into an ri−1 × ni ri matrix A for which
A AT = Iri−1

where ut is the (pM + 1) × 1 vector with entries
1, u1 (t), u2 (t), . . . , up (t), . . . , u1 (t − M + 1), u2 (t − M +
1), . . . , up (t − M + 1).

applies.
3

Tensor description of MIMO Volterra systems
The MIMO Volterra tensor consists of l (pM + 1)d entries,
which can quickly become practically infeasible to store as
M and d grow. We therefore propose to store the Volterra
tensor V in the TN format V (1) , . . . , V (d) , where the first
TN-core V (1) has dimensions l × pM + 1 × r1 . The other
TN-cores have sizes ri−1 × pM + 1 × ri , with rd = 1 for the
last core. The TN format reduces to a TT for the MISO case
(l = 1). This change in representation reduces the storage
requirement to O((d−1)(pM +1)r2 +(pM +1)lr). Lemma
2.1 then immediately tells us how to simulate the output
samples at time t as

To keep the notation simple, we first consider the following
discrete-time SISO Volterra system of degree d, namely,
y(t) = y0 (t) + y1 (t) + y2 (t) + · · · + yd (t),
= h0 +

d
X

M
−1
X

i=1 k1 ,...,ki =0

hi (k1 , . . . , ki )

i
Y

u(t − kj ),

j=1

(10)
where y(t), u(t) are the scalar output and input at time t respectively, M is the memory length and hi (k1 , . . . , ki ) denotes the ith Volterra kernel. Previous work that uses tensors in Volterra series describes each Volterra kernel as a
separate symmetric tensor [7]. Realizing from (10) that y(t)
is a multivariate polynomial in u(t), . . . , u(t − M + 1), we
define the SISO Volterra tensor as follows.

y(t) = (V (1) ×2 uTt ) (V (2) ×2 uTt ) · · · (V (d) ×2 uTt ),
(12)
with a computational complexity of O(d(pM +1)r+dr3 ). In
[1] a much faster simulation complexity of O(dRN log N )

4

is called a set of persistent exciting inputs of order d.

for computing N samples is obtained by using a symmetric polyadic representation for each of the Volterra kernels
separately. However, the method only works for the SISO
case and the identification of such a representation can be
problematic, since computation of the canonical rank R is
an NP-hard problem [12], which can be ill-posed [5].
4

We will from here on always assume that the inputs are
persistent exciting and that N ≥ rank(U ). The truncated
Volterra series for one particular output is a polynomial
in

+d
pM variables of degree d and therefore has pM
coeffipM
cients. A solution V (1) of the matrix equation (14) that con
+d
sists of only l pM
distinct entries therefore must correpM
spond with l symmetric tensors, where each column of V (1)
corresponds with one particular d-way symmetric tensor. As
a result, persistent exciting inputs and sufficient number of
samples N together with Lemma 4.1 has the following consequences:

MIMO Volterra system identification

In what follows we always consider the MIMO case. Observe now that we can rewrite (11) as
y(t)T = (ut d )T (V (1) )T ,

(13)

(1) the identification problem (14) has an infinite number
of solutions,
(2) the l unique minimal norm solutions of (14) correspond
with l symmetric tensors,
(3) the unicity of the symmetric solutions implies that no
other symmetric solutions exist,
(4) any solution V (1) of (14) can be turned into l minimal
norm solutions by symmetrizing each of the tensors
given by the columns of V (1) .

where V (1) is the Volterra tensor reshaped into a l × (pM +
1)d matrix. Writing out (13) for t = 0, 1, . . . , N − 1 leads
to the following matrix equation
Y = U (V (1) )T ,

(14)

where

T
Y := y(0) y(1) y(2) · · · y(N − 1) ,

T
U := u0 d u1 d u2 d · · · uN −1 d
.

Ideally, one would solve (14) for the minimal norm solutions, which acts as a regularization of the problem such that
the solution is uniquely defined. The two iterative methods
which we describe in the next subsections do not guarantee convergence to the minimal norm solutions. However, in
practice we observe that the norms of the obtained solutions
are quite close to the minimal ones.

Essentially, MIMO Volterra system identification is in its
most basic form solving the matrix equation (14). This immediately reveals the difficulty, as the size of the U matrix
is N × (pM + 1)d , which quickly becomes prohibitive even
for SISO systems and moderate degree d. This motivates us
to solve the following problem.

4.1

Problem 1 For a given set of measured time series
y1 , . . . , yl , u1 , . . . , up , memory M and degree d, solve the
matrix equation (14) for V in the TN format.

Alternating Linear Scheme method

The first method that we derive is the Alternating Linear
Scheme (ALS) method. The key idea of this method is to fix
the TN-ranks r1 , r2 , . . . , rd−1 and choose a particular initial
guess for V (1) , . . . , V (d) . Each of the TN-cores is then updated separately in an iterative fashion until convergence has
been reached. Once the core V (d) is updated, the algorithm
“sweeps” back towards V (1) and so on. It turns out that updating one of the cores is equivalent with solving a much
smaller linear system, which can be done very efficiently.
Before describing the form of the reduced linear system, we
first introduce the following notation

Before presenting the two numerical algorithms that solve
Problem 1, we first give an upper bound on the rank of U ,
together with some important implications.
Lemma
 4.1 The rank of the matrix U is upper bounded by
pM +d
pM .
Proof 1 Each row of U consists of (pM +1)d entries formed
d
by the repeated
 Kronecker product uk . There are however
pM +d
only pM distinct entries.

vk−1 := (V (1) ×2 uTt ) · · · (V (k−1) ×2 uTt ) ∈ Rl×rk−1 ,
vk+1 := (V (k+1) ×2 uTt ) · · · (V (d) ×2 uTt ) ∈ Rrk ×1 .

Lemma 4.1 motivates us to define a particular set of inputs
such that this upper bound is achieved.
Theorem 4.1 For outputs y1 , . . . , yl described by a MIMO
Volterra system (12) we have that

Definition 4.1 A set of input signals u1 , u2 , . . . , up such
that


pM + d
rank(U ) =
pM

T
y(t) = (vk+1
⊗ uTt ⊗ vk−1 ) vec(V (k) ).

5

(15)

Proof 2 We first rewrite (12) as

downloaded from https://github.com/kbatseli/
MVMALS.

y(t) = vk−1 (V (k) ×2 uTt ) vk+1 .

Algorithm 1 MIMO Volterra ALS Identification
Input: N samples of inputs u1 , . . . , up , outputs y1 , . . . , yl ,
degree d, memory M , TN ranks r1 , . . . , rd−1
Output: TN-cores V (1) , . . . , V (d) that solve Problem 1
Initialize right orthogonal TN-cores V (1) , . . . , V (d) of
prescribed ranks
while termination criterion not satisfied do
for i = 1, . . . , d − 1 do
vec(V (i) ) ← Compute and solve (16)
Vi ← reshape(V (i) , [ri−1 , (pM + 1)ri ])
Compute thin QR decomposition of Vi
V (i) ← reshape(Q, [ri−1 , (pM + 1), ri ])
Vi+1 ← reshape(V (i+1) , ri , (pM + 1)ri+1 ])
V (i+1) ← reshape(RVi+1 , [ri , (pM + 1), ri+1 ])
end for
Repeat the above loop in the reverse order
end while

This equation holds since it is the product of the l × rk−1
matrix vk−1 with the rk−1 × rk matrix (V (k) ×2 uTt ) with
the rk × 1 vector vk+1 , resulting in the l × 1 vector y(t).
We then have that
y(t) = vk−1 (V (k) ×2 uTt ) vk+1 ,
T
= (vk+1
⊗ vk−1 ) vec(V (k) ×2 uTt ),
T
= (vk+1
⊗ vk−1 ) vec(V (k) ×1 Irk−1 ×2 uTt ×3 Irk ),
T
= (vk+1
⊗ vk−1 ) (Irk ⊗ uTt ⊗ Irk−1 ) vec(V (k) ),
T
= (vk+1
⊗ 1 ⊗ vk−1 ) (Irk ⊗ uTt ⊗ Irk−1 ) vec(V (k) ),
T
= (vk+1
⊗ uTt ⊗ vk−1 ) vec(V (k) ),

where the second equation is obtained from (7), the fourth
equation is obtained from using (5) and the final equation
follows from (4). This concludes the proof.

A few remarks on the ALS method are in order.

The importance of Theorem 4.1 lies in the fact that it describes how the reduced linear system to update V (k) can
be constructed. Indeed, suppose we want to update V (k)
and keep all other TN-cores fixed. By using (15) for t =
0, . . . , N − 1 the following linear system is obtained
vec(Y T ) = Uk vec(V (k) ),

• The iterations can be terminated when the solution does
not exhibit any further improvement or when a certain
fixed maximal number of sweeps has been executed. In
our implementation we set a tolerance on the relative
residual ||Y − Ŷ ||2 /||Y ||2 , where Ŷ is the simulated
output from the obtained solution.
• It is proved in [11, 21] that under certain conditions the
ALS method enjoys strictly monotonous convergence.
Convergence to the unique minimal norm solutions cannot be guaranteed however.
• Also in [11, p. A701], it is proved that the QR decomposition step ensures that the condition number of each Uk
matrix in (16) is upper bounded by the condition number
of the large U matrix. This ensures the numerical stability
of the ALS method.

(16)

where Uk is a lN × rk−1 (pM + 1)rk matrix. Next to the
obvious savings in computational complexity compared to
solving (14), this effectively requires much fewer samples
N to perform the identification. Indeed, one only has to
make sure that lN ≥ rk−1 (pM + 1)rk . Computing the minimal norm solution of (16) requires the computation of the
pseudo-inverse of Uk , which requires O(4N l(rk−1 (pM +
1)rk )2 + 8(rk−1 (pM + 1)rk )3 ) computations.
Numerical stability of the ALS algorithm is guaranteed by
the use of an orthogonalization step [11, p. A690]. The key
idea is that all TN-cores are initialized to be right orthogonal
and are kept orthogonal during each step. After updating
V (1) by solving (16), this tensor is reshaped into a r0 (pM +
1) × r1 matrix from which a thin QR decomposition [9,
p. 230] is computed. This takes O(2l(pM + 1)r12 + 2r13 /2)
flops. The orthogonal matrix Q is then chosen as a new left
orthogonal V (1) TN-core. The corresponding R matrix is
then absorbed into the V (2) core by reshaping the core into a
r1 ×(pM +1)r2 matrix V2 and computing R V2 . Next, V (2)
is updated and orthogonalized, after which V (3) is updated
and so forth. When the iterative “sweep” has reached the end
of the TN, it reverses direction and in a similar way all cores
are made right orthogonal by QR decomposition. The whole
algorithm is presented as pseudocode in Algorithm 1. A
Matlab/Octave implementation of Algorithm 1 can be freely

While the ALS method has a small computational complexity, it suffers from the problem that all TN-ranks need to
be specified a priori. Finding a good choice can be quite
difficult when d is large. This is the main motivation for
the development of the Modified Alternating Linear Scheme
(MALS) algorithm, which not only updates the TN-cores
but also adapts the TN-ranks during each step. An additional
benefit is that the TN-cores are guaranteed to be either left
or right orthogonal so that no stabilizing QR step is required
anymore. These benefits, however, come at the cost of a
higher computational complexity.
4.2

Modified Alternating Linear Scheme method

The main idea of the MALS is in fact a simple one, namely,
to update two TN-cores V (k) , V (k+1) at a time by considering them as one “super-core” and keeping all other cores

6

where U , V are orthogonal matrices and S is a diagonal matrix with positive entries s1 ≥ . . . ≥ sq , with q =
min(rk−1 (pM +1), (pM +1)rk+1 ). Its computation requires
2
O(4(rk−1
(pM + 1)3 rk+1 ) + 8(rk−1 (pM + 1)3 rk+1 )2 +
3
9((pM + 1)3 rk+1
) flops. The numerical rank rk can be determined from a given tolerance τ such that s1 ≥ · · · ≥
srk ≥ τ ≥ · · · ≥ sq . If the sweep is going from left to right
then we update V (k) , V (k+1) as

fixed. The updated super-core is then decomposed into one
orthogonal and one non-orthogonal part, which are used as
updates for both V (k) and V (k+1) . This decomposition step
is achieved by computing a singular value decomposition
(SVD) and it is here where the TN-rank rk is updated. This
updating procedure is repeated in the same sweeping fashion
as with the ALS method until the solution has converged.
We now derive the reduced linear system for computing a
new super-core. The entries of the super-core V (k,k+1) are
defined as the contraction of V (k) with V (k+1) by summing
over all possible rk values
(k,k+1)

V i1 [i2 i3 ]i4 =

rk
X

(k)

V (k) := reshape(U , [rk−1 , (pM + 1), rk ]),
V (k+1) := reshape(SV T , [rk , (pM + 1), rk+1 ]),
which makes V (k) left orthogonal. If the sweep is going
from right to left then the updates are

(k+1)

V i1 i2 j V ji3 i4 .

j=1

V (k) := reshape(U S, [rk−1 , (pM + 1), rk ]),

The square brackets around i2 i3 indicate that these indices
need to be interpreted as one single multi-index such that
V (k,k+1) is an rk−1 × (pM + 1)2 × rk+1 3-way tensor.
Observe that the formation of the super-core has removed
all information on rk . It is now straightforward to verify that

V (k+1) := reshape(V T , [rk , (pM + 1), rk+1 ]),
where now V (k+1) is right orthogonal. The pseudocode
for the MALS algorithm is given in Algorithm 2. A Matlab/Octave implementation of Algorithm 2 can also be freely
downloaded from https://github.com/kbatseli/
MVMALS.

(V (k) ×2 uTt )(V (k+1) ×2 uTt ) = V (k,k+1) ×2 (uTt ) 2
holds. Using the same notation as in the ALS method we
now derive the reduced linear system for the MALS.

Algorithm 2 MIMO Volterra MALS Identification
Input: N samples of inputs u1 , . . . , up , outputs y1 , . . . , yl ,
degree d, memory M
Output: TN-cores V (1) , . . . , V (d) that solve Problem 1
Initialize right orthogonal TN-cores V (1) , . . . , V (d) of
ranks 1
while termination criterion not satisfied do
for i = 1, . . . , d − 1 do
vec(V (i,i+1) ) ← Compute and solve (18)
Vi,i+1 ← reshape(V (i,i+1) , [ri−1 (pM +1), (pM +
1)ri+1 ])
Compute the SVD of Vi,i+1
Determine numerical rank ri
V (i) ← reshape(U , [ri−1 , (pM + 1), ri ])
V (i+1) ← reshape(SV T , [ri , (pM + 1), ri+1 ])
end for
Repeat the above loop in the reverse order
end while

Theorem 4.2 For outputs y1 , . . . , yl described by a MIMO
Volterra system (12) we have that
T
y(t) = (vk+2
⊗ (uTt ) 2 ⊗ vk−1 ) vec(V (k,k+1) ).

(17)

Proof 3 By rewriting (12) this time as
y(t) = vk−1 (V (k,k+1) ×2 (uTt ) 2 ) vk+2 ,
it can be clearly seen that the rest of the proof is now identical to the ALS case.
Application of Theorem 4.2 for t = 0, . . . , N −1 now results
in the reduced linear system
vec(Y T ) = Uk,k+1 vec(V (k,k+1) ),

The same remarks as in the ALS case apply. The tolerance τ
for the determination of the numerical rank rk can be chosen
such that the error is below a certain threshold. In our implementation we opted for the default choice in Matlab/Octave,
which is τ =  s1 max(ri−1 (pM +1), (pM +1)ri+1 ), where
 is the machine precision. If the tolerance on the obtained
solution is set too low then the MALS algorithm will have
the tendency to increase the TN-ranks to very high values,
resulting in higher computational complexity. However, it
does not make sense in system identification to require that
the solution interpolates the measured output to a high accuracy (say up to the machine precision), which is essentially
overfitting.

(18)

where Uk,k+1 is an lN × rk−1 (pM + 1)2 rk+1 matrix.
The minimal norm solution of (18) can by computed in
O(4N l(rk−1 (pM + 1)2 rk+1 )2 + 8(rk−1 (pM + 1)2 rk+1 )3 )
flops.
Once the new super-core has been computed, it can be reshaped into an rk−1 (pM +1)×(pM +1)rk+1 matrix Vk,k+1 .
The SVD of Vk,k+1 is
Vk,k+1 = U S V T ,

7

5

Table 1
Run times, maximal TT-rank and number of estimated Volterra
tensor elements for an increasing degree d.
d
Run time [seconds]
max TT-rank (pM + 1)d

Numerical Experiments

In this section we demonstrate the two proposed identification algorithms. All computations were done on an Intel i5
quad-core processor running at 3.3 GHz with 16 GB RAM.
We are not aware of any other publicly available algorithms
that are able to identify high-degree MIMO Volterra systems, say, at d = 10.
5.1

Decaying multi-dimensional exponentials

First, we demonstrate the validity of Algorithms 1 and 2 by
means of an artificial SISO example. Symmetric Volterra
kernels were generated up to degree d = 10 for a fixed
memory M = 7 and containing exponentially decaying coefficients. The ith symmetric Volterra kernel hi contains the
entries

U †y

ALS

MALS

2

0.017

0.240

0.077

6

64

3

0.253

0.223

0.251

8

512

4

39.36

1.771

0.415

8

4096

5

NA

4.288

0.587

8

32768

6

NA

3.065

0.791

8

262144

7

NA

6.783

0.961

8

2097152

8

NA

13.11

1.199

8

16777216

9

NA

14.42

1.384

8

134217728

10

NA

18.37

1.576

8

1.0737e9

hi (k1 , . . . , ki ) = exp (−k12 − k22 − · · · − ki2 ).
Each of the ki indices attain the values 0, 0.1, 0.2, . . . , 0.6.
For each degree d a random input signal of 5000 samples, uniformly distributed over the interval [0, 1], was generated. The Volterra kernel was estimated by solving (14)
directly using the pseudoinverse of U where possible. The
MALS algorithm was used to determine a solution for which
||y − ŷ||2 /||y||2 < 1e−4 was satisfied. The TT-ranks determined by the MALS algorithm were then used to determine a solution with the same relative accuracy using
the ALS method. Both the MALS and ALS algorithms always ran with the first 700 samples of the input and output.
Table 1 lists the run times in seconds for the three methods, the maximal TT-rank and the total number of identified
Volterra kernel elements for each degree. Starting from degree d = 5 it was not possible anymore to obtain the Volterra
tensor directly since this would require the inversion of a
32768 × 32768 matrix. From d = 3, the maximal TT-rank
stabilizes to 8 for all TT-cores. Even though the ALS method
exhibits low computational complexity, its convergence is
much slower than MALS, resulting in larger run times. The
obtained solutions were validated by simulating the output
using the 4300 remaining input points. Fig. 3 shows both
the real and simulated output from the MALS solution for
samples 1640 up to 1700. The two output signals are almost indistinguishable. No difference between the ALS and
MALS solutions could be observed.
5.2

Fig. 3. Real and simulated output from the MALS solution with
d = 10, M = 7.

are added to the measured RF output, generating signals with
signal-to-noise (SNR) ratios ranging from 11dB up to 25dB.
The first 700 samples of the inputs and the noisy output are
then used to identify an M = 2, d = 11 two-input oneoutput Volterra system using the ALS method. The Volterra
kernel consists in this case of 9765625 entries. The TT-ranks
are all fixed to 2M + 1 = 5. The identified models were
then used to simulate the remaining 4300 samples of the
output. The SNR of the simulated output was computed by
comparing the simulated output with the original noiseless
output. Table 2 lists the SNR of the signals used in the
identification (ID SNR), the relative residual of the simulated
output, the run time of the identification in seconds and the
SNR of the simulated signal (SIM SNR). As expected, a
gradual improvement of the relative residual can be seen
as the SNR of the signals used for identification increases.
Although the residual remains high throughout the different
SNR levels, the SNR of the simulated output is much better,
with a consistent increase of 11dB. The run time varies
between 2 and 6 seconds. Fig. 4 shows the simulated output

Double balanced mixer

In this example we consider a double balanced mixer used
for upconversion. The output radio-frequency (RF) signal is
determined by a 100Hz sine low-frequency (LO) signal and
a 300Hz square-wave intermediate-frequency (IF) signal. A
phase difference of π/8 is present between the LO and IF
signals. All time series were sampled at 5 kHz for 1 second.
We investigate the effect of additive output noise on the
identified models. We define 5 different noise levels which

8

6

This article presented two new and remarkably efficient
identification algorithms for high-order MIMO Volterra systems. The identification problem was rephrased in terms of
the Volterra tensor, which is never explicitly constructed but
instead always stored in the highly economic TT format.
Both proposed identification algorithms are iterative, starting from an initial orthogonal guess for the TT-cores and
updating them until a desired accuracy is acquired. The algorithms are guaranteed to monotonically converge and numerical stability is ensured by retaining orthogonality in the
TT-cores. The efficiency of both identification algorithms
was demonstrated by numerical examples, where reliable
MIMO Volterra systems of degrees 10 were estimated in
only a few seconds. Even though its computational complexity is lower, the ALS method was found to converge slower
than the MALS algorithm when producing solutions with
the same accuracy. The MALS method showed a tendency
to increase the TT-ranks under the presence of high noise
levels. Extending the robustness of these methods to noisy
data will be the subject of further research.

Fig. 4. Real and simulated output from the ALS identified models
under different SNRs.

on the validation data for three Volterra models identified
under three different SNR levels (11dB, 16dB, 25 dB).
Next, the MALS algorithm was run on the same data to
also identify M = 2, d = 11 Volterra models. Based on
the ALS results we set the tolerance on the relative residual
to 0.5. This resulted in all TT-ranks being 5 for all cases.
Table 3 lists the SNR of the signals used in the identification
(ID SNR), the relative residual of the simulated output, the
run time of the identification in seconds and the SNR of
the simulated signal (SIM SNR). The relative residuals and
SNRs of the simulated output are very close to the results
obtained by the ALS identification. Just as in Example 5.1,
MALS is able to finish the identification faster than the ALS
method. Furthermore, the run time does not vary as much
as in the ALS case. Lowering the tolerance for the data with
low SNR resulted in the MALS method increasing the TT
ranks significantly, up to the point that the 5000 samples
were not sufficient anymore. This indicates a tendency of
the MALS method to overfit.
Table 2
ALS identification for 5 different SNR levels.
11dB 13dB 16dB 19dB
ID SNR

References
[1] K. Batselier, Z. Chen, H. Lui, and N. Wong. A tensor-based Volterra
series black-box nonlinear system identification and simulation
framework. In 2016 IEEE/ACM International Conference on
Computer-Aided Design (ICCAD), 2016.
[2] R. J. G. B. Campello, G. Favier, and W. C. do Amaral. Optimal
expansions of discrete-time Volterra models using Laguerre functions.
Automatica, 40(5):815 – 822, 2004.
[3] J. D. Carroll and J. J. Chang. Analysis of individual differences
in multidimensional scaling via an n-way generalization of “EckartYoung” decomposition. Psychometrika, 35(3):283–319, 1970.
[4] C.H. Cheng and E.J. Powers. Optimal Volterra kernel estimation
algorithms for a nonlinear communication system for PSK and QAM
inputs. IEEE Transactions on Signal Processing, 49(1):147–163,
2001.
[5] V. de Silva and L.H. Lim. Tensor rank and the ill-posedness of
the best low-rank approximation problem. SIAM J. on Matrix Anal.
Appl., 30(3):1084–1127, 2008.

25dB

||y−ŷ||2
||y||2

.255

.208

.151

.105

.052

Run time

2.3s

5.3s

6.4s

3.2s

2.2s

SIM SNR

22dB

24dB

27dB

30dB

37dB

[6] C. Diouf, M. Telescu, P. Cloastre, and N. Tanguy. On the Use
of Equality Constraints in the Identification of Volterra-Laguerre
Models. IEEE Signal Processing Letters, 19(12):857–860, Dec 2012.
[7] G. Favier, A. Y. Kibangou, and T. Bouilloc. Nonlinear system
modeling and identification using Volterra-PARAFAC models. Int.
J. Adapt. Control Signal Process, 26(1):30–53, January 2012.
[8] Doyle FJ III, Ronald K Pearson, and Babatunde A Ogunnaike.
Identification and control using Volterra models. Springer Science
& Business Media, 2012.

Table 3
MALS identification for 5 different SNR levels with tolerance of
0.5.
11dB 13dB 16dB 19dB 25dB
ID SNR
||y−ŷ||2
||y||2

.251

.214

.152

.105

.055

Run time

1.3s

1.1s

1.1s

1.1s

1.1s

SIM SNR

24dB

25dB

28dB

31dB

37dB

Conclusions

[9] G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns
Hopkins University Press, 3rd edition, October 1996.
[10] R. A. Harshman. Foundations of the PARAFAC procedure: Models
and conditions for an “explanatory” multi-modal factor analysis.
UCLA Working Papers in Phonetics, 16(1):84, 1970.
[11] S. Holtz, T. Rohwedder, and R. Schneider. The alternating linear
scheme for tensor optimization in the tensor train format. SIAM
Journal on Scientific Computing, 34(2):A683–A713, 2012.

9

[12] J. Håstad. Tensor rank is NP-complete. Journal of Algorithms,
11(4):644–654, 1990.
[13] Y. Kajikawa. Subband parallel cascade Volterra filter for linearization
of loudspeaker systems. In 2008 16th European Signal Processing
Conference, pages 1–5, Aug 2008.
[14] T. Katayama.
Subspace Methods for System Identification.
Communications and Control Engineering. Springer London, 2005.
[15] T. Kolda and B. Bader. Tensor decompositions and applications.
SIAM Review, 51(3):455–500, 2009.
[16] M.J. Korenberg and I.W. Hunter. The identification of nonlinear
biological systems: Volterra kernel approaches. Annals of biomedical
engineering, 24(2):250–268, 1996.
[17] E. Mumolo and D. Francescato. Adaptive predictive coding of
speech by means of Volterra predictors. In IEEE Winter Workshop
on Nonlinear Digital Signal Processing, 1993, pages 2.1.4.1–2.1.4.4,
1993.
[18] R. Orus. A Practical Introduction to Tensor Networks: Matrix Product
States and Projected Entangled Pair States. Annals Phys., 349:117–
158, 2014.
[19] I. Oseledets. Tensor-train decomposition. SIAM Journal on Scientific
Computing, 33(5):2295–2317, 2011.
[20] I. V. Oseledets. DMRG approach to fast linear algebra in the TT–
format. Comput. Meth. Appl. Math., 11(3):382–393, 2011.
[21] T. Rohwedder and A. Uschmajew. On local convergence of
alternating schemes for optimization of convex problems in the tensor
train format. SIAM Journal on Numerical Analysis, 51(2):1134–1162,
2013.
[22] L. Tan and J. Jiang. Adaptive Volterra filters for active control of
nonlinear noise processes. IEEE Transactions on Signal Processing,
49(8):1667–1676, Aug 2001.
[23] L. R. Tucker. Some mathematical notes on three-mode factor analysis.
Psychometrika, 31(3):279–311, 1966.
[24] P. Wambacq and W.M. Sansen. Distortion Analysis of Analog
Integrated Circuits. Kluwer Academic Publishers, Norwell, MA,
USA, 1998.
[25] S R. White. Density matrix formulation for quantum renormalization
groups. Phys. Rev. Lett., 69:2863–2866, Nov 1992.
[26] N. Wiener, J.A. Stratton, and M.I.O. Technology. Nonlinear Problems
in Random Theory. Literary Licensing, LLC, 2013.

10

