# A model reduction method for large-scale linear multidimensional dynamical
systems

M.A. Hamadi 22footnotemark: 2 Université Mohammed VI polytechnique, Lot 660,
Hay Moulay Rachid, Ben Guerir, 43150 Maroc<EMAIL_ADDRESS><EMAIL_ADDRESS>K. Jbilou11footnotemark: 1 Université du Littoral Côte
d’Opale, 50 Rue F. Buisson, BP 699–62228 Calais cedex, France; jbilou@univ-
littoral.fr A. Ratnani11footnotemark: 1

###### Abstract

In this work, we explore the application of multilinear algebra in reducing
the order of multidimentional linear time invariant (MLTI) systems. We use
tensor Krylov subspace methods as key tools, which involve approximating the
system solution within a low-dimensional subspace. We introduce the tensor
extended block and global Krylov subspaces and the corresponding Arnoldi based
processes. Using these methods, we develop a model reduction using projection
techniques. We also show how these methods could be used to solve large-scale
Lyapunov tensor equations that are needed in the balanced truncation method
which is a technique for order reduction. We demonstrate how to extract
approximate solutions via the Einstein product using the tensor extended block
Arnoldi and the extended global Arnoldi processes.

###### keywords:

Einstein product, Model recuction, Multi-linear dynamical systems, Tensor
Krylov subspaces.

In model reduction for high multi-linear time invariant (MLTI) systems,
multilinear algebra can be used to simplify the system representation and
reduce the computational complexity. This is achieved by decomposing the
system tensor into smaller, lower-dimensional tensors, which can be more
efficiently processed. One can use tensor decomposition techniques, such as
the Tucker, the higher-order singular value decomposition (HOSVD) or tensor
train decomposition (TTD) to decompose the system into simpler components,
which can then be used to identify and reduce the number of parameters needed
to describe the system.

Tensors are multi-dimensional arrays that can be used to represent higher-
order relationships between variables in a system, making them a powerful tool
for system analysis and reduction. The goal of model reduction in Multi-Linear
Time Invariant (MLTI) systems is to find a reduced-order representation that
accurately captures the behavior of the original system, while reducing the
computational cost. The resulting reduced-order model retains the essential
features of the original system and can be used for tasks such as simulation,
optimization and control. The use of multilinear algebra in this context can
significantly improve the efficiency and accuracy of the model reduction
process. In a various scientific and engineering applications, the dynamical
model over space and time can be described by a systems of PDEs, e.g., Navies-
Stokes equations in fluid dynamic, Schrodinguer equations in quantum mechanics
and heat equations in thermal transfer. Naturally, a spatial/temporal
discretization to such equations yields to tensors that evolve in time. A
straightforward idea that one can apply, is to transform tensors to matrices
or vectors, and then use some linear time invariant approaches as, projection
techniques onto Krylov subspaces [5, 14, 19, 25], Balanced Truncation (BT)
[16, 30, 31] or use some non-intrusive methods, such as the Loewner framework,
Eigensystem Realization Algorithm (ERA) [24] and Dynamic Mode Decomposition
(DMD) [36]. The effectiveness of those methods have been clearly proved in a
large body of literature. Other materials and details about model order
reduction techniques can be found in [2, 3, 7] and the references therein.
However the transformation from tensors to vectors leads to a computational
challenge due to high number of states/model parameters of the original system
that scales exponentially with the number of dimensions. The use of tensors
has become increasingly significant in many control applications, such as
state-space representation, system identification and optimal control, where
tensors can be used to represent the state of a control system in a compact
and efficient manner, making it easier to analyse the system and design
control strategies. The main idea is to use multilinear algebra to tensorize
the vector-based dynamic system representation into an equivalent tensor
representation. This strategy of guarding the compact tensor based
representations/computations have been established in many applications, such
as linear parameters varying models [4], nonlinear circuits via tensor
decomposition [29], solution of Lyapunov equations [32] and in simulation of
MLTI systems [27].

In this paper we are interested in projecting the MLTI system onto some Krylov
subspaces of a finite dimension using Einstein product. This technique of
projection helps to reduce the complexity of the original system while
preserving its essential features, and what we end up with is a lower-
dimensional reduced MLTI system that has the same structure as the original
one. The MLTI systems have been introduced in [11], where the states, outputs
and inputs are preserved in a tensor format and the dynamic is supposed to be
described by some multilinear operators. A technique named tensor unfolding
[9] which transforms tensor to matrix, allows to extend many concepts and
notions in the LTI theory to the tensor case. One of them is the use of the
transfer function that describes the input-output behaviour of the MLTI
system, and to quantify the efficiency of the projected reduced MLTI system.
Many properties in the theory of LTI systems such as, reachibility,
controllability and stability, have been generalized in the case of MLTI
systems based via the Einstein product [13]. A variety of applications in
model reduction and identification from tensor time series data have been
established in [34, 39]. The authors have shown the efficiency behind the
compact representation based on tensor algebra. A Balanced Truncation (BT)
model reduction for MLTI systems based on tensor structure has been introduced
in [13], where authors employed the tensor train decomposition (TTD) to make
computation more feasible. Through the paper, we will be interested in
projection and balanced truncation techniques using Krylov tensor based
subspaces, namely tensor extended global and block Krylov subspace methods
using new Arnoldi-based process. This will be done using the Einstein product
and multilinear algebra to decompose the system into simpler components via
projection. In the matrix case, the efficiency of the extended Krylov susbpace
method has been shown and proved in many fields, in LTI model reduction [1,
5], also for solving large scale algebraic matrix equations [22, 21, 37, 38].
We introduce tensor extended block and global Krylov subspaces and we use them
via an Arnoldi process in model reduction for MLTI, also for solving some
tensor equations based on Einstein product.

Recently, tensor Krylov subspace methods have gained interest beyond
scientific community. These methods are based on the idea of approximating the
solution of the system within a low-dimensional subspace, known as the Krylov
subspace, generated by applying a few tensor-tensor products. Iterative Krylov
subspace methods in tensor structure have shown efficiency and accurate
solution of high-dimensional problem. They can be used to solve large multi-
linear tensor equations [17, 18] with applications to color images and video
processing, using Einstein and T-products. These methods exploit the fact that
the computational cost grows linearly with the size of the problem. Another
application to solve multilinear systems based on Einstein product via tensor
inversion has been achieved in [9]. In a second part of the paper, we will be
interested in some large scale tensor equations such as the discrete Lyapunov
tensor equations that appears when using balanced truncation model order
reduction methods. We show how to extract approximate solution using tensor-
based Krylov subspace methods.

The paper is organized as follows. In Section 1, we present general notations
and give the main definitions used throughout the paper. An introduction to
the MLTI systems along with a description to the tensor classic and extended
global Krylov subspace based to construct high fidelity reduced MLTI systems
is shown in Section 2. The reduction process via projection using tensor
Krylov subspaces is established in Section 3. A description to the tensor
balanced truncation method is presented in Section 4. Computing approximate
solution to large scale discrete Lyapunov tensor equation via the tensor block
and global Krylov subspaces is described in Section 5. Finally, in Section 6,
some numerical experiments are assessed to show the accuracy of the proposed
methods.

##  1 Background and definitions

In this section, we provide with some notations and definitions that will be
used throughout this paper. Unless otherwise specified, we represent tensors
by calligraphic letters, we use lowercase and capital letters to identify
vectors and matrices, respectively. An $N$-mode tensor is a multidimensional
array and can be represented by
$\mathcal{X}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}}$, whose elements are
denoted by $x_{j_{1},\ldots,j_{N}}$ with $1\leq j_{i}\leq J_{i}$,
$i=1,\ldots,N$. If $N=0$, then $\mathcal{X}\in\mathbb{R}$ is a $0$th order
tensor that we call a scalar, a vector is a $1$th order tensor and a matrix is
$2$th order tensor. Further definitions and explanations about tensors could
be found in [26]. The tensor Einstein product is a multidimensional
generalizations of matrix product; for more details see [9, 13]

######  Definition 1.

Let ${\cal A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}}$ and in ${\cal
B}\in\mathbb{R}^{K_{1}\times\ldots\times K_{M}\times I_{1}\times\ldots\times
I_{L}}$ be two tensors. The Einstein product ${\cal A}\ast_{M}{\cal
B}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times I_{1}\times\ldots\times
I_{L}}$ is defined as follows

$({\cal A}\ast_{M}{\cal B})_{j_{1}\ldots j_{N}i_{1}\ldots
i_{L}}=\sum_{k_{1}=1}^{K_{1}}\sum_{k_{2}=1}^{K_{2}}\ldots\sum_{k_{M}=1}^{K_{M}}a_{j_{1}\ldots
j_{N}k_{1}\ldots k_{M}}b_{k_{1}\ldots k_{M}i_{1}\ldots i_{L}}.$ (1)

Some special tensor that will be considered in this paper : A nonzero tensor
${\cal D}=(d_{j_{1}\ldots j_{N}i_{1}\ldots
i_{N}})\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
J_{1}\times\ldots\times J_{N}}$, is called a Diagonal tensor if all the
entries $d_{j_{1}\ldots j_{N}i_{1}\ldots i_{N}}$ are equal to zero expect the
diagonal entries denoted by $d_{j_{1}\ldots j_{N}j_{1}\ldots j_{N}}\neq 0$. In
the case where all the diagonal entries are equal to $1$ (i.e.,
$d_{j_{1}\ldots j_{N}j_{1}\ldots j_{N}}=1$) then ${\cal D}$ is called the
identity tensor denoted by ${\cal I}$. Let ${\cal
A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times K_{1}\times\ldots\times
K_{M}}$ and ${\cal B}\in\mathbb{R}^{K_{1}\times\ldots\times K_{M}\times
J_{1}\times\ldots\times J_{N}}$ two tensors verified $b_{k_{1}\ldots
k_{M}j_{1}\ldots j_{N}}=a_{j_{1}\ldots j_{N}k_{1}\ldots k_{M}}$, then ${\cal
B}$ is called the transpose of ${\cal A}$ and denoted by ${\cal A}^{T}$. The
inverse of a square tensor is defined as follows.

######  Definition 2.

A square tensor ${\cal A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
J_{1}\times\ldots\times J_{N}}$ is invertible iff there exists a tensor
$\mathcal{X}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
J_{1}\times\ldots\times J_{N}}$ such that

${\cal A}\ast_{N}\mathcal{X}=\mathcal{X}\ast_{N}{\cal A}={\cal
I}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times J_{1}\times\ldots\times
J_{N}}.$

In that case, $\mathcal{X}$ is the inverse of ${\cal A}$ and is denoted by
${\cal A}^{-1}.$

######  Proposition 1.

Consider ${\cal A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}}$ and ${\cal
B}\in\small\mathbb{R}^{K_{1}\times\ldots\times K_{M}\times
J_{1}\times\ldots\times J_{N}}$. Then we have the following relations

  1. 1.

$({\cal A}\ast_{M}{\cal B})^{T}={\cal B}^{T}\ast_{N}{\cal A}^{T}.$

  2. 2.

${\cal I}_{M}\ast_{M}{\cal B}={\cal B}$ and ${\cal B}\ast_{N}{\cal
I}_{N}={\cal B}$.

The trace denoted by $tr(\cdot)$ of a square-order tensor ${\cal
A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times J_{1}\times\ldots\times
J_{N}}$ is given by

$tr({\cal A})=\sum_{j_{1}\ldots j_{N}}a_{j_{1}\ldots j_{N}j_{1}\ldots j_{N}}.$
(2)

The inner product of the two tensors
$\mathcal{X},\mathcal{Y}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}}$ is defined as

$<\mathcal{X},\mathcal{Y}>=tr(\mathcal{X}^{T}\ast_{N}\mathcal{Y}),$

where $\mathcal{Y}^{T}\in\mathbb{R}^{K_{1}\times\ldots\times K_{M}\times
J_{1}\times\ldots\times J_{N}}$ is the transpose tensor of $\mathcal{Y}$. If
$<\mathcal{X},\mathcal{Y}>=0$, we say that $\mathcal{X}$ and $\mathcal{Y}$ are
orthogonal. The corresponding norm is the tensor Frobenius norm given by

$\|\mathcal{X}\|=\sqrt{tr(\mathcal{X}^{T}\ast_{N}\mathcal{X})}.$ (3)

We recall the definition of the eigenvalues of a tensor mentioned in [28, 40].

######  Definition 3.

Let ${\cal A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
J_{1}\times\ldots\times J_{N}}$ be a tensor. The complex scalar $\lambda$
satisfied

${\cal A}\ast_{N}\mathcal{X}=\lambda\mathcal{X},$

is called an eigenvalue of ${\cal A}$, where
$\mathcal{X}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}}$ is a nonzero tensor
named an eigentensor. The set of all the eigenvalues of ${\cal A}$ is denoted
by $\Lambda({\cal A})$.

###  1.1 Tensor unfolding

Tensor unfolding plays an important role in tensor computations. It provides
simplification that allows to extended many notions to the tensor case that
are already presented in the matrix case. We reorganize the entries of the
tensor and present them in a vector or in a matrix format.

######  Definition 4.

Consider the following transformation $\Psi:\mathbb{T}_{j_{1}\ldots
j_{N}k_{1}\ldots k_{M}}\rightarrow\mathbb{M}_{|J|,|K|}$ with $\Psi({\cal
A})=A$ defined component wise as

${\cal A}_{j_{1}\ldots j_{N}k_{1}\ldots k_{M}}\rightarrow
A_{i\text{vec}(j,J),i\text{vec}(k,K)},$

where we refer to $J=\\{J_{1},\ldots,J_{N}\\}$, $K=\\{K_{1},\ldots,K_{M}\\}$
and $\mathbb{T}_{j_{1}\ldots j_{N}k_{1}\ldots k_{M}}$ is the set of all
tensors in $\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}}$. $\mathbb{M}_{|J|,|K|}$ is the set of matrices
in $\mathbb{R}^{|J|\times|K|}$ where $|J|=J_{1}\ldots J_{N}$ and
$|K|=K_{1}\ldots K_{N}$. The index mapping $i\text{vec}(\cdot,\cdot)$
introduced in [33] is given by

$i\text{vec}(j,J)=j_{1}+\sum_{l=2}^{N}(j_{l}-1)\prod_{l=1}^{s-1}J_{s}.$
$i\text{vec}(k,K)=k_{1}+\sum_{l=2}^{M}(k_{l}-1)\prod_{l=1}^{s-1}K_{s}.$

###  1.2 Block tensors

In this subsection, and analogously to block matrices, we describe briefly how
our tensors can be stored properly in a one tensor. To clarify our description
we need to start by giving the definition of an even-order paired tensor.

######  Definition 5.

A tensor is called even-order paired tensor if its order is always even, i.e.,
($2N$ for example), and the elements of the tensor are indicating using a
pairwise index notation, i.e., $a_{j_{1}i_{1}\ldots j_{N}i_{N}}$ for ${\cal
A}\in\mathbb{R}^{J_{1}\times I_{1}\times\ldots\times J_{N}\times I_{N}}$.

A definition of $n$-mode row block tensor, of two even-order paired tensor of
the same size is given by

######  Definition 6.

([13]) Consider ${\cal A},{\cal B}\in\mathbb{R}^{J_{1}\times
I_{1}\times\ldots\times J_{N}\times I_{N}}$ two even-order paired tensors, the
$n$-mode row block tensor denoted by $\displaystyle\big{|}{\cal A}\,\,{\cal
B}\big{|}_{n}\in\mathbb{R}^{J_{1}\times I_{1}\times\ldots\times J_{n}\times
2I_{n}\times\ldots\times J_{N}\times I_{N}}$ is defined by the following

$(\big{|}{\cal A}\,\,{\cal B}\big{|}_{n})_{j_{1}\ldots j_{N}i_{1}\ldots
i_{N}}=\left\\{\begin{array}[]{lll}a_{j_{1}\ldots j_{N}i_{1}\ldots
i_{N}},&&j_{k}=1,\ldots,J_{k},\,i_{k}=1,\ldots,I_{k}\,\forall k\\\
b_{j_{1}\ldots j_{N}i_{1}\ldots i_{N}},&&j_{k}=1,\ldots,J_{k},\,\forall
k,\,i_{k}=1,\ldots,I_{k}\,\,\text{for}\,\,k\neq n.\\\
&&\text{and}\,\,i_{k}=I_{k}+1,\ldots,2I_{k}\,\,\text{for}\,\,k=n.\end{array}\right.$
(4)

In this work, we are interested in the $4th$-order tensors (i.e., $N=2$), and
not especially paired ones (i.e., the paired wise index notation is not used).
We have noticed that the definition above stills valid for the case of
$4th$-order paired or non-paired tensor. In this regard, we specify more in
the following this block tensor notation by giving some examples. We suppose
now that ${\cal A}$ and ${\cal B}$ are two $4th$-order tensors in the space
$\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times I_{2}}$. By following the two
definitions above, we refer to $1$-mode row block tensor by $\big{|}{\cal
A}\,\,{\cal B}\big{|}_{1}\in\mathbb{R}^{J_{1}\times 2J_{2}\times I_{1}\times
I_{2}}$. We use the MATLAB colon operation $:$ to explain how to extract the
two tensors ${\cal A}$ and ${\cal B}$ as follows

$\big{|}{\cal A}\,\,{\cal B}\big{|}_{1}(:,1:J_{2},:,:)={\cal
A}\quad\text{and}\quad\big{|}{\cal A}\,\,{\cal
B}\big{|}_{1}(:,J_{2}+1:2J_{2},:,:)={\cal B}.$

The notation of $2$-mode row block tensor $\big{|}{\cal A}\,\,{\cal
B}\big{|}_{2}\in\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times 2I_{2}}$
described below is the mostly used in this paper. The extraction goes as

$\big{|}{\cal A}\,\,{\cal B}\big{|}_{2}(:,:,:,1:K_{2})={\cal
A}\quad\text{and}\quad\big{|}{\cal A}\,\,{\cal
B}\big{|}_{2}(:,:,:,K_{2}+1:2K_{2})={\cal B}.$

######  Remark 1.

We notice that if we consider $N=1$ (i.e., ${\cal A}$ and ${\cal B}$ are now
matrices, we refer to them as $A$ and $B$), then the notation $\big{|}{\cal
A}\,\,{\cal B}\big{|}_{1}=\big{|}{\cal A}\,\,{\cal
B}\big{|}_{2}=[A,B]\in\mathbb{R}^{J_{1}\times 2I_{1}}$ is now defined by the
standard block matrices definition.

Similarly to the $1$ or $2$-mode row block tensor, we can define the $1$ or
$2$-mode column block tensor for the $4th$-order tensors based on the
definition proposed by Can et al., in [11]. Using Matlab notation,
$\begin{vmatrix}{\cal A}\\\ {\cal B}\end{vmatrix}_{1}$ is the 1-mode column
block tensor in $\mathbb{R}^{2J_{1}\times J_{2}\times I_{1}\times I_{2}}$,
defined by

$\begin{vmatrix}{\cal A}\\\ {\cal B}\end{vmatrix}_{1}(1:J_{1},:,:,:)={\cal
A}\quad\text{and}\quad\begin{vmatrix}{\cal A}\\\ {\cal
B}\end{vmatrix}_{1}(J_{1}+1:2J_{1},:,:,:)={\cal B},$

and $\begin{vmatrix}{\cal A}\\\ {\cal B}\end{vmatrix}_{2}$ is 2-mode column
block tensor in $\mathbb{R}^{J_{1}\times J_{2}\times 2I_{1}\times I_{2}}$
given by

$\begin{vmatrix}{\cal A}\\\ {\cal B}\end{vmatrix}_{1}(:,:,1:I_{1},:)={\cal
A}\quad\text{and}\quad\begin{vmatrix}{\cal A}\\\ {\cal
B}\end{vmatrix}_{1}(:,:,I_{1}+1:2I_{1},:)={\cal B}.$

Based on the Einstein product, we can easily prove the following proposition
that will be mainly used during the computational process described in the
following sections.

######  Proposition 2.

Consider four tensors of the same size ${\cal A},{\cal B},{\cal C},{\cal
D}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}$, then we get

$\begin{array}[]{lll}\underbrace{\big{|}{\cal A}\,\,{\cal
B}\big{|}_{2}}_{\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times
2K_{2}}}\ast_{2}\underbrace{\big{|}{\cal C}\,\,{\cal
D}\big{|}_{1}}_{\in\mathbb{R}^{K_{1}\times 2K_{2}\times K_{1}\times
K_{2}}}&=&{\cal A}\ast_{2}{\cal C}+{\cal B}\ast_{2}{\cal
D}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}.\\\\[5.69046pt]
\underbrace{\big{|}{\cal A}\,\,{\cal
B}\big{|}_{1}}_{\in\mathbb{R}^{K_{1}\times 2K_{2}\times K_{1}\times
K_{2}}}\ast_{2}\underbrace{\big{|}{\cal C}\,\,{\cal
D}\big{|}_{2}}_{\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times
2K_{2}}}&=&\Big{|}\big{|}{\cal A}\ast_{2}{\cal C}\,\,\,{\cal A}\ast_{2}{\cal
D}\big{|}_{2},\,\,\big{|}{\cal B}\ast_{2}{\cal C}\,\,\,{\cal B}\ast_{2}{\cal
D}\big{|}_{2}\Big{|}_{1}\in\mathbb{R}^{K_{1}\times 2K_{2}\times K_{1}\times
2K_{2}}.\end{array}$

Let us consider that we have a $m$ tensors ${\cal
A}_{k}\in\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times I_{2}}$, we can apply
the definition associated to the $2$-mode row block tensor in succession to
create a one tensor denoted by ${\cal A}$ in the space
$\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times mI_{2}}$ and has the
following form

${\cal A}=\big{|}{\cal A}_{1}\quad{\cal A}_{2}\quad\ldots\quad{\cal
A}_{m}\big{|}\in\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times mI_{2}}.$ (5)

For more general definition we refer the reader to [[13], Definition $4.2$].
The process to create such tensor is defined as follows

  * -

If $m$ is even then we consider the following

$\mathcal{X}_{12}=\underbrace{\big{|}{\cal A}_{1}\quad{\cal
A}_{2}\big{|}_{2}}_{\in\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times
2I_{2}}},\quad\mathcal{X}_{34}=\underbrace{\big{|}{\cal A}_{3}\quad{\cal
A}_{4}\big{|}_{2}}_{\in\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times
2I_{2}}},\quad\ldots,\quad\mathcal{X}_{m-1m}=\underbrace{\big{|}{\cal
A}_{m-1}\quad{\cal A}_{m}\big{|}_{2}}_{\in\mathbb{R}^{J_{1}\times J_{2}\times
I_{1}\times 2I_{2}}},$

then we stored them two by two as follows

$\mathcal{X}_{1234}=\underbrace{\big{|}\mathcal{X}_{12}\quad\mathcal{X}_{34}\big{|}_{2}}_{\in\mathbb{R}^{J_{1}\times
J_{2}\times I_{1}\times
4I_{2}}},\quad\ldots,\quad\mathcal{X}_{m-3m-2m-1m}=\underbrace{\big{|}\mathcal{X}_{m-3m-2}\quad\mathcal{X}_{m-1m}\big{|}_{2}}_{\in\mathbb{R}^{J_{1}\times
J_{2}\times I_{1}\times 4I_{2}}},$

We keep repeating the process until the last 2-mode row block tensor obtained
and simply we presented in the following notation

${\cal A}=\big{|}{\cal A}_{1}\quad{\cal A}_{2}\quad\ldots\quad{\cal
A}_{m}\big{|}\in\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times mI_{2}}.$ (6)

  * -

If $m$ is odd, then we first choose $m-1$ tensors and applied the process
described above. We denote by $\mathcal{X}$ the obtained tensor in
$\mathbb{R}^{J_{1}\times J_{2}\times I_{1}\times(m-1)I_{2}}$ and then we
obtain the desired tensor ${\cal A}_{m}$ given by

${\cal A}=\big{|}\mathcal{X}\quad{\cal A}_{m}\big{|}\in\mathbb{R}^{J_{1}\times
J_{2}\times I_{1}\times mI_{2}}.$

An explanation of the MATLAB colon $:$ for the constructed tensor ${\cal A}$
goes as follows

${\cal A}_{l}={\cal
A}(:,:,:,(l-1)I_{2}+1:lI_{2})\quad\text{for}\,\,l=1,\ldots,m.$

In Sections 2 and 5, we will deal with some tensors in
$\mathbb{R}^{K_{1}\times mK_{2}\times K_{1}\times mK_{2}}$. Let us denote this
tensor $\mathcal{M}_{m}$, and suppose that $m=3$ (i.e.,
$\mathcal{M}_{3}\in\mathbb{R}^{K_{1}\times 3K_{2}\times K_{1}\times 3K_{2}}$),
by following the definitions above, the tensor $\mathcal{M}_{3}$ can be
presented as follows

$\mathcal{M}_{3}=\bigg{|}\,\Big{|}\,\big{|}\mathcal{M}_{1,1}\quad\mathcal{M}_{1,2}\big{|}_{2}\quad\mathcal{M}_{1,3}\Big{|}_{2},\,\,\Big{|}\,\big{|}\mathcal{M}_{2,1}\quad\mathcal{M}_{2,2}\big{|}_{2}\quad\mathcal{M}_{2,3}\Big{|}_{2},\,\Big{|}\,\big{|}\mathcal{M}_{3,1}\quad\mathcal{M}_{3,2}\Big{|}_{2}\quad\mathcal{M}_{3,3}\Big{|}_{2}\,\,\bigg{|}_{1},$
(7)

where $\mathcal{M}_{i,j}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times
K_{2}}$ and the three $2$-mode row block tensor given above have been
constructed as mentioned in (5), where we concatenate every tensor
$\mathcal{M}_{i,j}$ successively. To facilitate the notation, we will refer to
$\mathcal{M}_{m}$ as

$\mathcal{M}_{3}=\begin{vmatrix}\mathcal{M}_{1,2}&\mathcal{M}_{1,2}&\mathcal{M}_{1,3}\\\
\mathcal{M}_{2,1}&\mathcal{M}_{2,2}&\mathcal{M}_{2,3}\\\
\mathcal{M}_{3,1}&\mathcal{M}_{3,2}&\mathcal{M}_{3,3}\end{vmatrix}\in\mathbb{R}^{K_{1}\times
3K_{2}\times K_{1}\times 3K_{2}}.$ (8)

The tensors $\mathcal{M}_{i,j}$ could be defined from the tensor
$\mathcal{M}_{m}$ by using the MATLAB colon as follows

$\mathcal{M}_{i,j}=\mathcal{M}_{3}(:,(i-1)K_{2}+1:iK_{2},:,(j-1)K_{2}+1:jK_{2})\quad\text{for}\,\,i,j=1,2,3.$

This could be followed by a generalisation to $m\in\mathbb{N}$ (i.e.,
$\mathcal{M}_{m}\in\mathbb{R}^{K_{1}\times mK_{2}\times K_{1}\times mK_{2}}$)
defined as follows

$\mathcal{M}_{m}=\begin{vmatrix}\mathcal{M}_{1,1}&\mathcal{M}_{1,2}&\ldots&\mathcal{M}_{1,m}\\\
\mathcal{M}_{2,1}&\mathcal{M}_{2,2}&\ldots&\mathcal{M}_{2,m}\\\
\vdots&\ddots&\ddots&\vdots\\\
\mathcal{M}_{m,1}&\mathcal{M}_{m,2}&\ldots&\mathcal{M}_{m,m}\end{vmatrix}\in\mathbb{R}^{K_{1}\times
mK_{2}\times K_{1}\times mK_{2}},$ (9)

and in the same way, we can describe every tensor $\mathcal{M}_{i,j}$ using
the MATLAB colon by the following

$\mathcal{M}_{i,j}=\mathcal{M}_{m}(:,(i-1)K_{2}+1:iK_{2},:,(j-1)K_{2}+1:jK_{2})\in\mathbb{R}^{K_{1}\times
K_{2}\times K_{1}\times K_{2}}\quad\text{for}\,\,i,j=1,2,\ldots,m.$

For a simple presentation, we suppose that $m=2$ and ${\cal A},{\cal B}$ are
two tensors in $\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}$, then
using the definitions described above, we can easy prove the following
proposition

######  Proposition 3.

  * 1.

$\begin{vmatrix}\mathcal{M}_{1,1}&\mathcal{M}_{1,2}\\\
\mathcal{M}_{2,1}&\mathcal{M}_{2,2}\end{vmatrix}=\Big{|}\,\big{|}\mathcal{M}_{1,1}\quad\mathcal{M}_{1,2}\big{|}_{2}\,\,\big{|}\mathcal{M}_{2,1}\quad\mathcal{M}_{2,2}\big{|}_{2}\Big{|}_{1}=\Big{|}\,\big{|}\mathcal{M}_{1,1}\quad\mathcal{M}_{2,1}\big{|}_{1}\,\,\big{|}\mathcal{M}_{1,2}\quad\mathcal{M}_{2,2}\big{|}_{1}\Big{|}_{2}$.

  * 2.

$\big{|}{\cal A}\quad{\cal
B}\big{|}_{2}\ast_{2}\begin{vmatrix}\mathcal{M}_{1,1}&\mathcal{M}_{1,2}\\\
\mathcal{M}_{2,1}&\mathcal{M}_{2,2}\end{vmatrix}=\big{|}{\cal
A}\ast_{2}\mathcal{M}_{1,1}+{\cal B}\ast_{2}\mathcal{M}_{2,1}\quad{\cal
A}\ast_{2}\mathcal{M}_{1,2}+{\cal B}\ast_{2}\mathcal{M}_{2,2}\big{|}_{2}.$

  * 3.

$\begin{vmatrix}\mathcal{M}_{1,1}&\mathcal{M}_{1,2}\\\
\mathcal{M}_{2,1}&\mathcal{M}_{2,2}\end{vmatrix}\ast_{2}\big{|}{\cal
A}\quad{\cal B}\big{|}_{1}=\big{|}\mathcal{M}_{1,1}\ast_{2}{\cal
A}+\mathcal{M}_{1,2}\ast_{2}{\cal B}\quad\mathcal{M}_{2,1}\ast_{2}{\cal
A}+\mathcal{M}_{2,2}\ast_{2}{\cal B}\big{|}_{1}.$

##  2 Multidimensional linear model reduction

Consider the following discrete time MLTI system

$\left\\{\begin{split}\mathcal{X}_{k+1}&={\cal A}\ast_{N}\mathcal{X}_{k}+{\cal
B}\ast_{M}{\cal U}_{k},\quad\mathcal{X}_{0}=0,\\\ \mathcal{Y}_{k}&={\cal
C}\ast_{N}\mathcal{X}_{k},\end{split}\right.$ (10)

in the continuous case, a continuous MLTI system could be described as follows

$\left\\{\begin{split}\dot{\mathcal{X}}(t)&={\cal
A}\ast_{N}\mathcal{X}(t)+{\cal B}\ast_{M}{\cal
U}(t),\quad\mathcal{X}_{0}=0,\\\ \mathcal{Y}(t)&={\cal
C}\ast_{N}\mathcal{X}(t),\end{split}\right.$ (11)

where $\dot{\mathcal{X}}(t)$ is the derivative of the tensor
$\mathcal{X}(t)\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}}$, ${\cal
A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times J_{1}\times\ldots\times
J_{N}}$ is a square tensor, ${\cal B}\in\mathbb{R}^{J_{1}\times\ldots\times
J_{N}\times K_{1}\times\ldots\times K_{M}}$ and ${\cal
C}\in\mathbb{R}^{I_{1}\times\ldots\times I_{M}\times J_{1}\times\ldots\times
J_{N}}$. The tensors ${\cal U}_{k}\in\mathbb{R}^{K_{1}\times\ldots\times
K_{M}}$ (${\cal U}(t)$ in the continuous case) and
$\mathcal{Y}_{k}\in\mathbb{R}^{I_{1}\times\ldots\times I_{M}}$
($\mathcal{Y}(t)$ in the continuous case) are the control and output tensors,
respectively.

######  Proposition 4.

The transfer function associated to the dynamical system (11) is given by

${\cal F}(s)={\cal C}\ast_{N}(s{\cal I}-{\cal A})^{-1}\ast_{N}{\cal B}.$ (12)

###### Proof.

Given the transformation $\Psi$ defined in 4, we can show that the multi-
linear system (10) could be transformed into an equivalent matrix linear
dynamical system, defined by

$\left\\{\begin{split}x_{k+1}&=Ax_{k}+Bu_{k},\quad x_{0}=0\\\
y_{k}&=Cx_{k},\end{split}\right.$ (13)

where $A,B$ and $C$ are the matrices resulting from the matricization of the
tensors ${\cal A},{\cal B}$ and ${\cal C}$ respectively by using $\Psi$. The
same remark goes for the vectors $x_{k},u_{k}$ and $y_{k}$. System of the form
(13), can be represented in a frequency domain via Z-transform by the
following transfer function (cf. [10])

$F(s)=C(sI-A)^{-1}B.$

From [[9], Lemma $3.1$], it is proved that for two matrices $X$ and $Y$ with
appropriate dimensions, we have $\Psi^{-1}(XY)=\mathcal{X}\ast\mathcal{Y}$.
Then, we get the following result

$\begin{split}{\cal F}(s)=\Psi^{-1}(F(s))&=\Psi^{-1}(C(sI-A)^{-1}B)\\\
&=\Psi^{-1}(C)\ast_{N}\Psi^{-1}((sI-A)^{-1})\ast_{N}\Psi^{-1}(B)\\\ &={\cal
C}\ast_{N}(s{\cal I}-{\cal A})^{-1}\ast_{N}{\cal B}.\end{split}$

∎

It is worthy to mention that the theoretical concepts that characterize the
linear dynamical system (13) such as stability, reachability and observability
have been extended to the MLTI systems. Originally, these theoretic concepts
have been introduced by Can et al., in [13] in the case of a discrete MLTI
system (10). The results shown below are exactly the ones described for the
matrix case if one considers the isomorphism $\Psi$ described previously. If
we consider a non-zero initial tensor state $\mathcal{X}_{0}$, and analogously
to the matrix case, the solution to the discrete MLTI systems (10) is given by

$\begin{array}[]{lll}\mathcal{X}_{k}&=&{\cal A}^{\star
k}\ast_{N}\mathcal{X}_{0}+\displaystyle\sum_{j=0}^{k-1}{\cal A}^{\star
k-j-1}\ast_{N}{\cal B}\ast_{M}{\cal U}_{j},\end{array}$ (14)

where ${\cal A}^{\star k}=\underbrace{{\cal A}\ast\ldots\ast{\cal A}}_{k}$.
The notion of stability of an MLTI system is given as follows

######  Definition 7.

Given the MLTI system (10). The system is called

  * 1.

Asymptotically stable if $\|\mathcal{X}_{k}\|\leq\kappa\|\mathcal{X}_{0}\|$
for some $\kappa>0$.

  * 2.

Stable if $\|\mathcal{X}_{k}\|\rightarrow 0$ as $t\rightarrow\infty$.

We have the following results related to the stability of an MLTI system.

######  Proposition 5.

([11])

  * 1.

System (10) is considered stable if and only if the magnitudes of all the
eigenvalues of ${\cal A}$ (Definition 3) are less than or equal to $1$.
Additionally, for any eigenvalues equal to $1$, the algebraic and geometry
multiplicity must be equal.

  * 2.

If the magnitudes of all the eigenvalues of ${\cal A}$ are less than $1$ than
System (10) is called asymptotically stable.

In the discrete case, the definitions of the reachability and the
observability of a discrete MLTI system are defined below.

######  Definition 8.

Consider the discrete MLTI system (10). The pair (${\cal A},{\cal B}$) is
reachable on $\left[k_{0},k_{1}\right]$ if and only if the reachability
Gramian

$\mathcal{P}_{d}(k_{0},k_{1})=\sum_{k=k_{0}}^{k_{1}-1}{\cal
A}^{k_{1}-k-1}\ast_{N}{\cal B}\ast_{M}{\cal B}^{T}\ast_{N}({\cal
A}^{k_{1}-k-1})^{T},$ (15)

is a weakly symmetric positive-definite square tensor, see [11] for more
details about this notion.

The reachability Gramian is a solution of the following discrete Lyapunov
(called also Stein) tensor equation

$\mathcal{P}_{d}-{\cal A}\ast_{N}\mathcal{P}_{d}\ast_{N}{\cal A}^{T}={\cal
B}\ast_{N}{\cal B}^{T}.$ (16)

Likewise, the observability of an MLTI system is given as follows.

######  Definition 9.

The pair (${\cal A},{\cal C}$) is observable on $\left[k_{0},k_{1}\right]$ if
and only if the observability Gramian

${\cal Q}_{d}(k_{0},k_{1})=\sum_{k=k_{0}}^{k_{1}-1}({\cal
A}^{k_{1}-k-1})^{T}\ast_{N}{\cal C}^{T}\ast_{M}{\cal C}\ast_{N}{\cal
A}^{k_{1}-k-1},$ (17)

is a weakly symmetric positive-definite square tensor, see [11] for more
details about this notion.

This observability Gramian ${\cal Q}_{d}$ determines how much each state
influences future outputs [35]. It fulfils the following discrete Discrete-
Lyapunov tensor equation

${\cal Q}_{d}-{\cal A}^{T}\ast_{N}{\cal Q}_{d}\ast_{N}{\cal A}={\cal
C}^{T}\ast_{M}{\cal C}.$ (18)

The two discrete Lyapunov tensor equations (16) and (18) play a crucial role
in the Balanced Truncation model order reduction for MLTI systems described in
[12]. A technique based on TT-algebra [12] and outer product has been
established to solve such equations. In the next section, we describe a tensor
Krylov subspace projection technique for solving properly such equations. As
in the case of LTI systems, in order to construct a reduced model using a
Krylov subspace techniques, the whole process remains on approximating the
associated transfer function. This could also be shared with the MLTI systems.
In the following, we introduce a tensor global Arnoldi Krylov method using
Einstein product to construct such approximation that yields at the end to the
reduced MLTI model. We use the classic Krylov subspace and also introduce the
tensor extended Krylov subspace.

###  2.1 Tensor based global Krylov subspace methods

####  2.1.1 Classic Krylov subspace processing

The first step in the approximation of the associated transfer function (12)
to the MLTI system (10) is that we rewrite it as

${\cal F}(s)={\cal C}\ast_{N}\mathcal{X},$

where $\mathcal{X}$ verify the following multi-linear system

$(s{\cal I}-{\cal A})\ast_{N}\mathcal{X}={\cal B}.$ (19)

In order to find an approximation to ${\cal F}(s)$, it remains to find an
approximation to the above multi-linear system, which can be done by using a
projection Krylov subpace technique. Let us first define the $m$-th tensor
Krylov subspace

${\cal K}_{m}({\cal A},{\cal B})=\text{span}\\{{\cal B},{\cal A}\ast{\cal
B},\ldots,{\cal A}^{\star m-1}\ast{\cal
B}\\}\subseteq\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}},$ (20)

where ${\cal A}^{\star k}=\underbrace{{\cal A}\ast\ldots\ast{\cal A}}_{k}$. We
restrict our attention in the case where $N=M=2$ for the sake of simplicity.
In other words, $4th$-order tensors will be taken into consideration in the
following sections. The results can be easily extended to the cases where
$N\geq 2$ and $M\geq 2$. In the sequel, unless otherwise specified, we always
refer to the Einstein product between two $4th$-order tensors $"\ast_{2}"$ as
$"\ast"$.
The tensor global Anroldi process is an analogue to the one described for
matrices [23]. The algorithm linked to it is outlined below, and it takes a
form that differs from the one described in [17].

Algorithm 1 Tensor global Arnoldi algorithm

  1. 1.

Input : ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$,
${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$ and a
fixed integer $m$.

  2. 2.

Output : $\widetilde{\mathcal{V}}_{m+1}\in\mathbb{R}^{J_{1}\times J_{2}\times
K_{1}\times(m+1)\,K_{2}}$ and $\overline{H}_{m}\in\mathbb{R}^{m+1\times m}$.

  3. 3.

Set $\beta=\|{\cal B}\|_{F}$ and $\mathcal{V}_{1}={\cal B}/\beta$. (MATLAB
notation $\widetilde{\mathcal{V}}_{m+1}(:,:,:,1:K_{2})=\mathcal{V}_{1}$)

  4. 4.

for $j=1\cdots m$

  5. 5.

$\mathcal{W}={\cal A}\ast\mathcal{V}_{j}$

  6. 6.

for $i=1,\ldots,j$.

  7. 7.

$h_{i,j}=<\mathcal{V}_{i},\mathcal{W}>$.

  8. 8.

$\mathcal{W}=\mathcal{W}-h_{i,j}\mathcal{V}_{i}$.

  9. 9.

end for.

  10. 10.

$h_{j+1,j}=\|\mathcal{W}\|_{F}$. If $h_{j+1,j}=0$, stop. Otherwise
$\mathcal{V}_{j+1}=\mathcal{W}/h_{j+1}.$

  11. (MATLAB notation. $\widetilde{\mathcal{V}}_{m+1}(:,:,:,jK_{2}+1:(j+1)K_{2})=\mathcal{V}_{j+1}$).

  12. 11.

end for.

What we end up with, after $m$ steps, is the following decomposition

${\cal
A}\ast\widetilde{\mathcal{V}}_{m}=\widetilde{\mathcal{V}}_{m+1}\ast(\overline{H}_{m}\circledast{\cal
I}_{K}),$ (21)

where $\widetilde{\mathcal{V}}_{m+1}\in\mathbb{R}^{J_{1}\times J_{2}\times
K_{1}\times(m+1)K_{2}}$ contains $\mathcal{V}_{i}\in\mathbb{R}^{J_{1}\times
J_{2}\times K_{1}\times K_{2}}$ for $i=1,\ldots,m+1$, as it is explained in
the algorithm above. $\overline{H}_{m}\in\mathbb{R}^{m+1\times m}$ is a
Hessenberg matrix of the following form

$\overline{H}_{m}=\begin{bmatrix}h_{1,1}&h_{1,2}&\ldots&h_{1,m}\\\
h_{2,1}&h_{2,2}&\ldots&h_{2,m}\\\ 0&h_{3,2}&\ldots&h_{3,m}\\\
\vdots&\ddots&\ddots&\vdots\\\ 0&\ddots&&h_{m,m}\\\
0&\ldots&0&h_{m+1,m}\end{bmatrix}.$

${\cal I}_{K}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}$ is the
identity tensor, a diagonal tensor with entries

${\cal
I}_{K}(i,j,k,l)=\delta_{i,k}\delta_{j,l},\quad\text{where}\quad\delta_{p,q}=1,\,\,\text{if}\,\,p=q,\,\,\delta_{p,q}=0,\,\,\text{otherwise}.$
(22)

The product $\circledast$ used in the decompositions is defined as follows.

######  Definition 10.

Given a matrix $P\in\mathbb{R}^{m\times n}$ and a tensor ${\cal
J}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}$, the resulted
tensor $\mathcal{R}=P\circledast{\cal J}$ is defined as follows

  1. 1.

$\mathcal{R}\in\mathbb{R}^{K_{1}\times mK_{2}\times K1\times nK_{2}}$.

  2. 2.

$\mathcal{R}(i,:,i,:)=P\otimes{\cal J}(i,:,i,:)$, for $i=1,\ldots,K_{1}$,
where $\otimes$ is the Kronecker product.

Notice that if ${\cal J}$ is a matrix in $\mathbb{R}^{K_{1}\times K_{2}}$ then
the product $\circledast$ becomes the Kronecker product $\otimes$. The
following result is used in the computational process described in the next
sections. It is summarized in the following straightforward proposition.

######  Proposition 6.

Let $\widetilde{\mathcal{V}}_{m}\in\mathbb{R}^{J_{1}\times J_{2}\times
K_{1}\times mK_{2}}$ be the tensor generated from Algorithm 1 Let $P,Q^{T}$
two matrices in $\mathbb{R}^{m\times n}$, ${\cal
I}_{K}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}$ the identity
tensor. Then we have

$\widetilde{\mathcal{V}}_{m}\ast(PQ\circledast{\cal
I}_{K})=\widetilde{\mathcal{V}}_{m}\ast((P\circledast{\cal
I}_{K})\ast(Q\circledast{\cal I}_{K})),$

Notice that the computed $\mathcal{V}_{i}$’s from Algorithm 1 form an
orthonormal basis of the tensor Krylov subspace 20; i.e.
$\widetilde{\mathcal{V}}_{m}^{T}\diamond\widetilde{\mathcal{V}}_{m}=I_{m}$,
where $I_{m}\in\mathbb{R}^{m\times m}$ is the identity matrix. The product
$\diamond$ is an analogous to the one defined in the global Arnoldi process
described for the matrix case in [23] as follows

$\widetilde{\mathcal{V}}_{m}^{T}\diamond\widetilde{\mathcal{V}}_{m}=\begin{bmatrix}<\mathcal{V}_{1},\mathcal{V}_{1}>&<\mathcal{V}_{1},\mathcal{V}_{2}>&\ldots&<\mathcal{V}_{1},\mathcal{V}_{m}>\\\
<\mathcal{V}_{2},\mathcal{V}_{1}>&<\mathcal{V}_{2},\mathcal{V}_{2}>&\ldots&<\mathcal{V}_{2},\mathcal{V}_{m}>\\\
\vdots&\ddots&\ddots&\ldots\\\
<\mathcal{V}_{m},\mathcal{V}_{1}>&\ldots&\ldots&<\mathcal{V}_{m},\mathcal{V}_{m}>\end{bmatrix}\in\mathbb{R}^{m\times
m},$ (23)

where $<\cdot,\cdot>$ is the inner-product defined in the previous section for
tensors. Thanks to [8] which concerns some properties of the $\diamond$
product in the matrix case, we can obtain similar results in the tensor case

######  Proposition 7.

Given tensors ${\cal E},{\cal F},\mathcal{G}\in\mathbb{R}^{J_{1}\times
J_{2}\times K_{1}\times mK_{2}}$ and matrix $S\in\mathbb{R}^{m\times m}$. Then
we have

  1. 1.

$({\cal E}+{\cal F})^{T}\diamond\mathcal{G}={\cal
E}^{T}\diamond\mathcal{G}+{\cal F}^{T}\diamond\mathcal{G}.$

  2. 2.

${\cal E}^{T}\diamond({\cal F}+\mathcal{G})={\cal E}^{T}\diamond{\cal F}+{\cal
E}^{T}\diamond\mathcal{G}$.

  3. 3.

${\cal E}^{T}\diamond({\cal F}\ast(S\circledast I_{K}))=({\cal
E}^{T}\diamond{\cal F})S$.

Based on these results, and using the definition of the product $\circledast$,
we can easily prove the following proposition that will be of a great use in
terms of simplification as we will show in the following sections. Let
$\widetilde{\mathcal{V}}_{m+1}$ be the orthonormal tensor given by Algorithm
1, then as for the matrix case, we have the classical algebraic relations

$\widetilde{\mathcal{V}}_{m+1}^{T}\diamond({\cal
A}\ast\widetilde{\mathcal{V}}_{m})=\overline{H}_{m},\quad\text{and}\quad\widetilde{\mathcal{V}}_{m}^{T}\diamond({\cal
A}\ast\widetilde{\mathcal{V}}_{m})=H_{m}.$ (24)

where $H_{m}\in\mathbb{R}^{m\times m}$ is also a Hessenberg matrix extracted
from $\overline{H}_{m}$ by deleting the last row.

####  2.1.2 The extended Krylov subspace process

For the matrix case, extended block and global Arnoldi methods are more
effective than the classical ones. This is usually the case In model order
reduction techniques [1, 5], for large scale dynamical systems; [22, 21, 37,
38]. In this section, we describe the extended tensor global Arnoldi process.
Before that let us give some definitions and remarks about the extended tensor
Krylov subspace. Let ${\cal A}$ and $\mathcal{V}$ be two tensors of
appropriate dimensions, then the $2m$-dimension extended tensor global Krylov
subspace denoted by ${\cal K}^{e}_{m}({\cal A},\mathcal{V})$, can be defined
as follows

$\begin{array}[]{lll}{\cal K}^{e}_{m}({\cal A},\mathcal{V})&=&{\cal
K}_{m}({\cal A},\mathcal{V})+{\cal K}_{m}({\cal A}^{-1},{\cal
A}^{-1}\ast\mathcal{V})\\\ &=&\text{span}\left\\{{\cal A}^{-\star
m}\ast\mathcal{V},\ldots,{\cal A}^{-1}\ast\mathcal{V},\mathcal{V},\ldots,{\cal
A}^{\star\,m-1}\ast\mathcal{V}\right\\},\end{array}$ (25)

where ${\cal A}^{\star k}=\underbrace{{\cal A}\ast\ldots\ast{\cal A}}_{k}$,
${\cal A}^{\star-k}=\underbrace{{\cal A}^{-1}\ast\ldots\ast{\cal A}^{-1}}_{k}$
and ${\cal K}_{m}({\cal A},\mathcal{V})$ is the tensor classic global Krylov
subspace defined in (20) associated to the pair $({\cal A},\mathcal{V})$. It
is worth noting that if we consider the isomorphism $\Psi$ defined previously
to define the matrices $\Psi({\cal A})=A\in\mathbb{R}^{J_{1}J_{2}\times
J_{1}J_{2}}$ and $\Psi(\mathcal{V})=V\in\mathbb{R}^{J_{1}J_{2}\times
K_{1}K_{2}}$, then we can generalize all the results obtained in the matrix
case, from the extended global Krylov subspace process applied to the pair
$(A,V)$, to the tensor structure via the Einstein product. Further details are
given below. Next, we describe the process for the construction of the tensor
$\widetilde{\mathcal{V}}_{2m}$ associated to the extended global Krylov
subspace defined above. The process can be guaranteed via the following
extended Arnoldi algorithm. As mentioned earlier, our interest is focused only
on $4th$-order tensors, but the following results remain true also for higher
order.

Algorithm 2 Tensor Extended Global Arnoldi (TEGA) algorithm

  1. 1.

Input : ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$,
${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$ and a
fixed integer $m$.

  2. 2.

Output : $\widetilde{\mathcal{V}}_{2(m+1)}\in\mathbb{R}^{J_{1}\times
J_{2}\times K_{1}\times 2(m+1)K_{2}}$ and
$\overline{H}_{m}\in\mathbb{R}^{2(m+1)\times 2m}$.

  3. 3.

Compute the global QR-decomposition of $|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times 2K_{2}}$ i.e.,
$\text{gQR}(|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2})=[\mathcal{V}_{1},\omega]$. (MATLAB notation
$\widetilde{\mathcal{V}}_{2(m+1)}(:,:,:,1:2K_{2})=\mathcal{V}_{1}$,
$\omega=\begin{bmatrix}\omega_{1,1}&\omega_{1,2}\\\
0&\omega_{2,2}\end{bmatrix}\in\mathbb{R}^{2\times 2}$)

  4. 4.

for $j=1\cdots m$

  5. 5.

Set $\mathcal{W}=|{\cal A}\ast\mathcal{V}_{j}^{1},{\cal
A}^{-1}\ast\mathcal{V}_{j}^{2}|_{2}$, where
$\left\\{\begin{array}[]{lll}\mathcal{V}_{j}^{1}&=&\mathcal{V}_{j}(:,:,:,(j-1)K_{2}+1:jK_{2}),\\\
\mathcal{V}_{j}^{2}&=&\mathcal{V}_{j}(:,:,:,jK_{2}+1:(j+1)K_{2}),\end{array}\right.$

  6. 6.

for $i=1,\ldots,j$.

  7. 7.

$H_{i,j}=\mathcal{V}_{i}^{T}\diamond\mathcal{W}\in\mathbb{R}^{2\times 2}$.

  8. 8.

$\mathcal{W}=\mathcal{W}-\mathcal{V}_{i}\ast(H_{i,j}\circledast{\cal I}_{K})$.

  9. 9.

end for.

  10. 10.

Compute the global QR (gQR) of $\mathcal{W}$, i.e.,
$\text{gQR}(\mathcal{W})=[\mathcal{V}_{j+1},H_{j+1,j}]$

  11. 11.

end for.

Algorithm 2 computes an orthonormal tensor basis corresponding to
$\widetilde{\mathcal{V}}_{2(m+1)}\in\mathbb{R}^{J_{1}\times J_{2}\times
K_{1}\times 2(m+1)K_{2}}$, associated to the extended global Krylov subspace
(25) applied to the pair $({\cal A},{\cal B})$. This means that for
$i,j=1,\ldots,m$, the following is verified

$\mathcal{V}_{i}^{T}\diamond\mathcal{V}_{j}=0_{2\times 2}\quad i\neq
j\quad\text{and}\quad\mathcal{V}_{i}^{T}\diamond\mathcal{V}_{i}=I_{2\times
2}.$ (26)

An upper Hessenber matrix $\overline{H}_{m}\in\mathbb{R}^{2(m+1)\times 2m}$ is
also generated with upper block matrices $H_{i,j}\in\mathbb{R}^{2\times 2}$.
In Step $3$ and Step $10$ of Algorithm 2, we compute a global QR decomposition
to get the element $\mathcal{V}_{j}$ of the tensor basis
$\widetilde{\mathcal{V}}_{2(m+1)}$. This decomposition is analogue to the one
described in [1, 5] for the extended global Krylov subspace in the matrix
case. We use this decomposition to get $|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2}=\mathcal{V}_{1}\ast(\omega\circledast{\cal I}_{K})$ in Step $3$ and
$\mathcal{W}=\mathcal{V}_{j+1}\ast(H_{j+1,j}\circledast{\cal I}_{K})$ in Step
$10$, where product $\circledast$ is the one described in Definition 10. From
steps $5$ to $9$, we obtain the following relations

$\begin{array}[]{lll}\mathcal{V}_{j+1}\ast(H_{j+1,j}\circledast
I_{K})&=&|{\cal A}\ast\mathcal{V}_{j}^{1},{\cal
A}^{-1}\ast\mathcal{V}_{j}^{2}|_{2}-\displaystyle\sum_{i=1}^{j}\mathcal{V}_{i}\ast(T_{i,j}\circledast{\cal
I}_{K}),\\\ H_{j+1,j}&=&\mathcal{V}_{j+1}^{T}\diamond|{\cal
A}\ast\mathcal{V}_{j}^{1},{\cal A}^{-1}\ast\mathcal{V}_{j}^{2}|_{2}\\\
&=&\mathcal{V}_{j+1}^{T}\diamond\displaystyle\sum_{i=1}^{j+1}\mathcal{V}_{i}\ast(H_{i,j}\circledast{\cal
I}_{K}).\end{array}$ (27)

After $m$ steps of Algorithm 2 we can show that

$\begin{array}[]{lll}{\cal
A}\ast\widetilde{\mathcal{V}}_{2m}&=&\widetilde{\mathcal{V}}_{2(m+1)}\ast(\overline{T}_{m}E_{m}^{T}\circledast{\cal
I}_{K})\\\ &=&\widetilde{\mathcal{V}}_{2m}\ast(T_{m}\circledast{\cal
I}_{k})+\mathcal{V}_{m+1}\ast(T_{m+1,m}E_{m}^{T}\circledast{\cal
I}_{K}),\end{array}$ (28)

where $E_{m}$ is the last $2m\times 2$ block column of the identity matrix
$I_{m}\in\mathbb{R}^{2m\times 2m}$,
$\overline{T}_{m}=\widetilde{\mathcal{V}}_{2(m+1)}^{T}\diamond({\cal
A}\ast\widetilde{\mathcal{V}}_{2m})\in\mathbb{R}^{2(m+1)\times 2m}$ and
$T_{m}$ is a $2m\times 2m$ block Hessenberg matrix defined by
$T_{m}=\widetilde{\mathcal{V}}_{2m}^{T}\diamond({\cal
A}\ast\widetilde{\mathcal{V}}_{2m})\in\mathbb{R}^{2m\times 2m}$, where
$T_{i,j}=\mathcal{V}_{i}^{T}\diamond({\cal A}\ast\mathcal{V}_{j})$, for
$i,j=1,\ldots,m$. We summarize in the following result a recursion to compute
$T_{m}$ by avoiding additional tensor products with ${\cal A}$.

######  Proposition 8.

([20]) Let
$\overline{H}_{m}=[h_{:,1},\ldots,h_{:,m}]\in\mathbb{R}^{2(m+1)\times 2m}$ be
the upper Hessenberg matrix generated from Algorithm 2 where
$h_{:,i}\in\mathbb{R}^{2(m+1)\times 1}$ is the $ith$-column. A simple
computation of the upper Hessenberg matrix
$\overline{T}_{m}=[t_{:,1},\ldots,t_{:,m}]\in\mathbb{R}^{2(m+1)\times 2m}$ is
as follows.
The odd columns of the $\overline{T}_{m}$ are such that

$t_{:,2j-1}=h_{:,2j-1},\quad\text{for}\,\,j=1,\ldots,m,$

and the even column are described as

  * -

$t_{:,2}=\frac{1}{\omega_{22}}(\omega_{11}e_{1}^{2(m+1)}-\omega_{12}h_{:,1}),$

  * -

$t_{:,2j+2}=\frac{1}{h_{2j+2,2j}}(e_{2j}^{2(m+1)}-t_{:,1:2j+1}h_{1:2j+1,2j})\quad\text{for}\,\,j=1,\ldots,m-1$,

where $e_{i}^{(k)}$ is the $ith$ column vector of the identity matrix $I_{k}$
and $\omega_{1,1},\omega_{1,2}$ and $\omega_{2,2}$ are as defined in step $2$
from Algorithm 2.

##  3 Order reduction via projection

Based on the tensor classic and extended Krylov subspace processes described
above, we are now able to provide a reduction process to get a reduced MLTI
system to the original one (10). We recall our original MLTI system

$\left\\{\begin{split}\mathcal{X}_{k+1}&={\cal A}\ast\mathcal{X}_{k}+{\cal
B}\ast{\cal U}_{k},\quad\mathcal{X}_{0}=0,\\\ \mathcal{Y}_{k}&={\cal
C}\ast\mathcal{X}_{k},\end{split}\right.$ (29)

where the tensors ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times
J_{2}}$ and ${\cal B},{\cal C}^{T}\in\mathbb{R}^{J_{1}\times J_{2}\times
K_{1}\times K_{2}}$. For simplicity, we consider product $\ast$ as the
Einstein product $\ast_{2}$ between $4th$-order tensors. The associated
transfer function mentioned previously in Proposition 4 is given as

${\cal F}(s)={\cal C}\ast(s{\cal I}-{\cal A})^{-1}\ast{\cal B}.$ (30)

We seek to construct the following reduced MLTI system

$\left\\{\begin{split}\hat{\mathcal{X}}_{k+1}&=\hat{{\cal
A}}\ast\hat{\mathcal{X}}_{k}+\hat{{\cal B}}\ast\hat{{\cal
U}}_{k},\quad\mathcal{X}_{0}=0,\\\ \hat{\mathcal{Y}}_{k}&=\hat{{\cal
C}}\ast\hat{\mathcal{X}}_{k},\end{split}\right.$ (31)

by using Proposition 4, we can describe the associated transfer function as
follows

$\hat{{\cal F}}(s)=\hat{{\cal C}}\ast(s\hat{{\cal I}}-\hat{{\cal
A}})^{-1}\ast\hat{{\cal B}}.$

In order to measure the accuracy of the resulting reduced system, one can
compute the error-norm $\|{\cal F}-\hat{{\cal F}}\|$ with respect to a certain
norm. Similarly to reduction process described in for the matrix case, an
analogue process using a tensorial structure could be given. We first consider
one of the Krylov subspace defined in the previous section, the tensor classic
or extended Krylov subspaces, where we project the initial MLTI system. The
two subspaces are of $m$-dimension and $2m$-dimension respectively. We
consider tensor $\mathcal{V}_{m}$ as the corresponding basis tensor to one of
the tensor Krylov subspace. After approximating the full order state
$\mathcal{X}_{k}$ by $\mathcal{V}_{m}\ast\hat{\mathcal{X}}_{k}$, and
analogously to the matrix case, by using the Petrove-Galerkin condition
technique, we end up with the following

$\left\\{\begin{split}\mathcal{V}_{m}^{T}\ast(\mathcal{V}_{m}\ast\hat{\mathcal{X}}_{k+1}&-{\cal
A}\ast\mathcal{V}_{m}\ast\hat{\mathcal{X}}_{k}-{\cal B}\ast{\cal U}_{k})=0,\\\
\mathcal{Y}_{k}&={\cal
C}\ast\mathcal{V}_{m}\ast\hat{\mathcal{X}}_{k},\end{split}\right.$ (32)

finally, we obtain the desired reduced MLTI system (31), with the following
tensorial structure

$\hat{{\cal A}}=\mathcal{V}_{m}^{T}\ast({\cal
A}\ast\mathcal{V}_{m}),\quad\hat{{\cal B}}=\mathcal{V}_{m}^{T}\ast{\cal
B},\quad\hat{{\cal C}}={\cal C}\ast\mathcal{V}_{m}.$ (33)

In the following we give explicitly the expression of the tensorial structure
(33) by means of the two processes based on the two Krylov subspaces already
mentioned. To reach the desired goal, we need to find an approximation
$\hat{{\cal F}}$ to ${\cal F}$, and this comes back to find an approximation
$\mathcal{X}_{m}$ to $\mathcal{X}$ that verifies the multi-linear system (19)
in the tensor classic Krylov subspace, i.e., $\mathcal{X}_{m}\in{\cal
K}_{m}({\cal A},{\cal B})$, or in the tensor extended Krylov subspace i.e.,
$\mathcal{X}_{m}\in{\cal K}^{e}_{m}({\cal A},{\cal B})$. The process to find
such approximation based on a tensor classic Krylov subspace (20) is described
as follows.

$\mathcal{X}_{m}=\widetilde{\mathcal{V}}_{m}\ast(y_{m}\circledast{\cal
I}_{K})\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}},$ (34)
$\mathcal{R}_{m}={\cal B}-(s{\cal I}-{\cal
A})\ast\mathcal{X}_{m}\,\perp\,{\cal K}_{m}({\cal A},{\cal B}),$ (35)

where $y_{m}$ is a vector with $m$-components and
$\widetilde{\mathcal{V}}_{m}$ is the orthonormal basis associated to the
tensor Krylov subspace (20). We know from the decomposition (21) that

${\cal
A}\ast\widetilde{\mathcal{V}}_{m}=\widetilde{\mathcal{V}}_{m}\ast(H_{m}\circledast{\cal
I}_{K})+h_{m+1,m}\mathcal{V}_{m+1}\ast(e_{m}^{T}\circledast{\cal I}_{K}),$

and by using the fact ${\cal B}=\|{\cal B}\|\mathcal{V}_{1}$ and Proposition
7, the equality (35) could be simplified as follows

$\begin{split}\widetilde{\mathcal{V}}_{m}^{T}\diamond({\cal B}-(s{\cal
I}-{\cal A})\ast(\widetilde{\mathcal{V}}_{m}\ast(y_{m}\circledast{\cal
I}_{K}))&=0,\\\ \widetilde{\mathcal{V}}_{m}^{T}\diamond{\cal
B}&=\widetilde{\mathcal{V}}_{m}^{T}\diamond((s\widetilde{\mathcal{V}}_{m}-{\cal
A}\ast\widetilde{\mathcal{V}}_{m})\ast(y_{m}\circledast{\cal I}_{K})),\\\
\|{\cal B}\|e_{1}^{m}&=(sI_{m}-H_{m})y_{m},\\\
y_{m}&=(sI_{m}-H_{m})^{-1}\|{\cal B}\|e_{1}^{m}.\end{split}$

Now, we can have an explicit expression of the approximation of ${\cal F}(s)$
denoted by ${\cal F}_{m}(s)$ and described as

$\begin{split}{\cal F}(s)\approx\hat{{\cal F}}(s)={\cal
C}\ast\mathcal{X}_{m}&={\cal
C}\ast(\widetilde{\mathcal{V}}_{m}\ast(y_{m}\circledast{\cal I}_{K}))\\\
&={\cal C}\ast(\widetilde{\mathcal{V}}_{m}\ast(((sI_{m}-H_{m})^{-1}\|{\cal
B}\|e_{1}^{m})\circledast{\cal I}_{K}))\\\ &={\cal
C}\ast\widetilde{\mathcal{V}}_{m}\ast(s{\cal I}_{K}^{m}-(H_{m}\circledast{\cal
I}_{K}))^{-1}\ast(\|{\cal B}\|e_{1}^{m}\circledast{\cal I}_{K}),\end{split}$
(36)

where ${\cal I}_{K}^{m}\in\mathbb{R}^{K_{1}\times mK_{2}\times K_{1}\times
mK_{2}}$ is the identity tensor. Now, from the approximation $\hat{{\cal
F}}(s)$, we can conclude the reduced MLTI system

$\left\\{\begin{split}\hat{\mathcal{X}}_{k+1}&=\hat{{\cal
A}}\ast\hat{\mathcal{X}}_{k}+\hat{{\cal B}}\ast{\cal U}_{k},\\\
\hat{\mathcal{Y}}_{k}&=\hat{{\cal
C}}\ast\hat{\mathcal{X}}_{k},\end{split}\right.$ (37)

with the associated tensor structure

$\hat{{\cal A}}=H_{m}\circledast{\cal I}_{K},\quad\hat{{\cal B}}=\|{\cal
B}\|e_{1}^{m}\circledast{\cal I}_{K},\quad\hat{{\cal C}}={\cal
C}\ast\widetilde{\mathcal{V}}_{m}.$ (38)

In order to quantify how well the approximated transfer function $\hat{{\cal
F}}(s)$ is closer to the original one ${\cal F}(s)$, we propose in the
following a simplified expression to the error-norm $\|{\cal F}-\hat{{\cal
F}}\|$.

######  Proposition 9.

Let ${\cal F}(s)$ and $\hat{{\cal F}}(s)$ be the two transfer functions
associated to the original MLTI (10) and the reduced one (37) respectively.
Based on the previous results, we get

$\|{\cal F}(s)-\hat{{\cal F}}(s)\|\leq\|{\cal
C}\ast\widetilde{\mathcal{V}}_{m}\ast(s{\cal I}-{\cal
A})^{-1}\|\ast\|h_{m+1,m}((e_{m}^{T}(sI_{m}-H_{m})^{-1})\circledast{\cal
I}_{K})\ast{\cal B}_{m}\|.$ (39)

###### Proof.

The error ${\cal F}(s)-\hat{{\cal F}}(s)$ could be expressed as

$\begin{split}{\cal F}(s)-\hat{{\cal F}}(s)&={\cal C}\ast(s{\cal I}-{\cal
A})^{-1}\ast{\cal B}-{\cal C}_{m}\ast(s{\cal I}_{K}^{m}-{\cal
A}_{m})^{-1}\ast{\cal B}_{m}\\\ &={\cal C}\ast(s{\cal I}-{\cal
A})^{-1}\|\ast\left[{\cal B}-(s{\cal I}-{\cal
A})\ast\widetilde{\mathcal{V}}_{m}\ast(s{\cal I}_{K}^{m}-{\cal
A}_{m})^{-1}\ast{\cal B}_{m}\right],\end{split}$

by using the decomposition (21), we obtain

$\begin{split}{\cal F}(s)-\hat{{\cal F}}(s)&={\cal C}\ast(s{\cal I}-{\cal
A})^{-1}\ast\left[{\cal B}-\underbrace{\widetilde{\mathcal{V}}_{m}\ast
B_{m}}_{{\cal B}}+h_{m+1,m}\mathcal{V}_{m+1}\ast(e_{m}^{T}\circledast{\cal
I}_{K})\ast(s{\cal I}_{K}^{m}-{\cal A}_{m})^{-1}\ast{\cal B}_{m}\right]\\\
&={\cal C}\ast(s{\cal I}-{\cal
A})^{-1}\ast\left[h_{m+1,m}\mathcal{V}_{m+1}\ast((e_{m}^{T}(sI_{m}-H_{m})^{-1})\circledast{\cal
I}_{K})\ast{\cal B}_{m}\right],\end{split}$

this concludes the proof. ∎

If we consider the tensor extended global Krylov subspace ${\cal
K}^{e}_{m}({\cal A},{\cal B})$ (25), we can get the approximated transfer
function $\hat{{\cal F}}$ by following the same process described above based
on the tensor classic global Krylov subspace, then, the approximated
$\hat{{\cal F}}$ transfer function has the following form

$\hat{{\cal F}}(s)=\hat{{\cal C}}\ast(s{\cal I}_{2m}-\hat{{\cal
A}})^{-1}\ast\hat{{\cal B}},$ (40)

with the associated structure matrices

$\hat{{\cal A}}=T_{m}\circledast{\cal I}_{K},\quad\hat{{\cal
B}}=(\widetilde{\mathcal{V}}_{2m}^{T}\diamond{\cal B})\circledast{\cal
I}_{K},\quad\hat{{\cal C}}={\cal C}\ast\widetilde{\mathcal{V}}_{2m}.$ (41)

From Algorithm 2, we know that

  * •

$|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2}=\mathcal{V}_{1}\ast(\omega\circledast{\cal I}_{K})$ where
$\omega\in\mathbb{R}^{2\times 2}$ is an upper triangular matrix from Algorithm
2, then,

${\cal B}_{m}=(\omega_{1,1}e_{1}^{(2m)})\circledast{\cal I}_{K},$

where $e_{1}^{(2m)}$ is the first column of the identity matrix $I_{2m}$.

  * •

$T_{m}=\widetilde{\mathcal{V}}_{2m}^{T}\diamond({\cal
A}\ast\widetilde{\mathcal{V}}_{2m})$ a $2m\times 2m$ is a Hessenberg matrix
that can be computed in a efficient way using Proposition 8 and
$\widetilde{\mathcal{V}}_{2m}$ is the tensor computed by Algorithm 2.

######  Remark 2.

It is worth noting that in the case of a continuous-time MLTI system described
as

$\left\\{\begin{split}\dot{\mathcal{X}}(t)&={\cal A}\ast\mathcal{X}(t)+{\cal
B}\ast{\cal U}(t),\quad\mathcal{X}_{0}=0,\\\ \mathcal{Y}(t)&={\cal
C}\ast\mathcal{X}(t),\end{split}\right.$ (42)

the process described above for constructing a discrete-time reduced MLTI
system can be also applied to the continuous-time MLTI system, we end up by
the following reduced system

$\left\\{\begin{split}\dot{\mathcal{X}_{m}}(t)&={\cal
A}_{m}\ast\mathcal{X}_{m}(t)+{\cal B}_{m}\ast{\cal U}_{t},\\\
\mathcal{Y}_{m}(t)&={\cal C}_{m}\ast\mathcal{X}_{m}(t),\end{split}\right.$
(43)

where the coefficient tensorial are as follows

${\cal A}_{m}=\mathcal{V}_{m}^{T}\ast({\cal A}\ast\mathcal{V}_{m}),\quad{\cal
B}_{m}=\mathcal{V}_{m}^{T}\ast{\cal B},\quad{\cal C}_{m}={\cal
C}\ast\mathcal{V}_{m}.$ (44)

##  4 Tensor Balanced Truncation

Besides Krylov subspace model reduction methods for large scale LTI systems,
another class of methods is the known Balanced Truncation method (BT)
introduced by Moore [31]. The main challenge in the BT is to solve large-scale
discrete-time Lyapunov equations in order to obtain the system Gramians that
will be used to generate a reduced model. In this section, we provide with a
tensor Balanced Truncation method. The process of this method is as follows.
First, we need to solve two discrete Lyapunov equations

$\mathcal{P}-{\cal A}\ast\mathcal{P}\ast{\cal A}^{T}={\cal B}\ast{\cal
B}^{T},\qquad{\cal Q}-{\cal A}^{T}\ast{\cal Q}\ast{\cal A}={\cal
C}^{T}\ast{\cal C},$ (45)

where ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$ and
${\cal B},{\cal C}^{T}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times
K_{2}}$. As mentioned before, $\mathcal{P}$ and ${\cal Q}$ are known as the
reachability and the observability Gramians. It is mentioned previously, that
$\mathcal{P}$ and ${\cal Q}$ are weakly symmetric positive-definite square
tensor, then we can obtain the Cholesky-like factors of the two gramians
described as follows

$\mathcal{P}={\cal Z}\ast{\cal Z}^{T},\qquad{\cal
Q}=\mathcal{W}\ast\mathcal{W}^{T}.$ (46)

Where the tensors ${\cal Z}$ and $\mathcal{W}$ are of appropriate dimension
presented in low-rank form. The next step consists in computing the singular
value decomposition (SVD) of $\Psi(\mathcal{W}^{T}\ast{\cal Z})=U\Sigma
V^{T}$, where $\Psi$ is the isomorphism defined in (4). A reasonable
computation is required behind this step since tensor $\mathcal{W}$ and ${\cal
Z}$ are not of large dimensions. The matrix $\Sigma$ is a diagonal matrix that
contains the singular values known also as the Hankel singular values of the
associated MLTI system [12]. This decomposition could be expressed as follows

$\begin{split}\Psi(\mathcal{W}^{T}\ast{\cal Z})&=U\Sigma V^{T}\\\
&=\begin{bmatrix}U^{1}\,\,U^{2}\end{bmatrix}\begin{bmatrix}\Sigma^{1}&\\\
&\Sigma^{2}\end{bmatrix}\begin{bmatrix}(V^{1})^{T}\\\
(V^{2})^{T}\end{bmatrix},\end{split}$ (47)

where $\Sigma_{1}\in\mathbb{R}^{r^{2}\times r^{2}}$ and
$\Sigma_{2}\in\mathbb{R}^{(N^{2}-r^{2})\times(N^{2}-r^{2})}$. As described in
[16, 30, 31], a truncation step could be established, and by truncating the
states that correspond to small Hankel singular values in $\Sigma_{2}$. The
tensorial structure of the reduced MLTI system is given as follows

${\cal A}_{r}=\mathcal{Y}_{r}^{T}\ast({\cal A}\ast\mathcal{X}_{r}),\quad{\cal
B}_{r}=\mathcal{Y}_{r}^{T}\ast{\cal B},\quad{\cal C}_{r}={\cal
C}\ast\mathcal{X}_{r},$ (48)

with

$\mathcal{Y}_{r}=\mathcal{W}_{2}\ast\Psi^{-1}(U^{1})\ast\Psi^{-1}(\Sigma_{1})^{(-1/2)},\qquad\mathcal{X}_{r}={\cal
Z}_{1}\ast\Psi^{-1}(V^{1})\ast\Psi^{-1}(\Sigma_{1})^{(-1/2)},$ (49)

where $\Psi^{(-1)}$ is the inverse of the isomorphism defined in (4).
Regarding the solutions $\mathcal{P}$ and ${\cal Q}$ of the two Lyapunov
equations (45), we suggest a solution in a factored form via the tensor
extended block and global Krylov subspace projection methods. These two
methods are the focus of the upcoming section.

##  5 Tensor disctete-time Lyapunov equations

In this section we discuss the process to get approximate solutions to the
tensor discrete Lyapunov equations of the form

$\mathcal{X}-{\cal A}\ast\mathcal{X}\ast{\cal A}^{T}={\cal B}\ast{\cal
B}^{T},$ (50)

where ${\cal A},{\cal B}$ are tensors of appropriate dimensions and
$\mathcal{X}$ is the unknown tensor. The solution of such equation is the main
task in Balanced Truncation model order reduction method for MLTI systems
described in [12]. They also play a crucial role in control theory, for
instance, the equation (50) is what we get after a discretization of the $2D$
heat equations with control [12, 32]. It is obvious that, if the dimension of
(50) is small, then one can transform (50) to a matrix Lyapunov equation and
using some efficient direct techniques described in [6, 21, 38]. But in the
case of large scale tensors, choosing a process based only on tensors will be
of benefits than using some matricization techniques that will cost us both
computationally and in terms of memory. In the next section, we propose a
process to solve discrete Lyapunov tensor equations by means of tensor Krylov
subspace techniques. We use the global process described previously which is
based on a tensor classic or extended Krylov subspace. Moreover, a block
process by means of the tensor extended block Krylov subspace for solving (50)
is also presented. A comparison between all the methods is assessed in the
numerical section.

###  5.1 Tensor global Krylov method for discrete-Lyapunov equation

We recall the tensor Lyapunov equation that we are interested in

$\mathcal{X}-{\cal A}\ast\mathcal{X}\ast{\cal A}^{T}={\cal B}\ast{\cal
B}^{T},$ (51)

where ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$,
${\cal B}\in\mathbb{R}^{K_{1}\times K_{2}\times K_{1}\times K_{2}}$ and
$\mathcal{X}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$ is the
unknown tensor. Analogously to the matrix discrete-Lyapunov equation [6], we
state that (51) has a unique solution if the following is satisfied

$\lambda_{j_{1}j_{2}j_{1}j_{2}}\,\lambda_{j_{1}j_{2}j_{1}j_{2}}\neq 1,\quad
j_{1}=1,\ldots,J_{1};\,\,j_{2}=1,\ldots,J_{2}$

where $\lambda_{j_{1}j_{2}j_{1}j_{2}}$ are the eigenvalues of ${\cal A}$.

####  5.1.1 The extended process

In the following we describe the case of constructing the approximation by
using a projection onto the tensor extended global Krylov subspace defined in
(25). The approximate solution is given by

$\mathcal{X}_{m}=\widetilde{\mathcal{V}}_{2m}\ast(Y_{2m}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T},$ (52)

where $\widetilde{\mathcal{V}}_{2m}$ is a tensor obtained after running
Algorithm 2 to the pair (${\cal A},{\cal B}$), ${\cal
I}_{2m}\in\mathbb{R}^{K_{2}\times 2mK_{2}\times K_{1}\times 2mK_{2}}$ is the
identity tensor and the product $\circledast$ is as defined in Definition 10.
$Y_{2m}\in\mathbb{R}^{2m\times 2m}$ is a small matrix obtained via the
following Galerkin orthogonality condition

$\widetilde{\mathcal{V}}_{2m}^{T}\diamond(\mathcal{R}_{m}\diamond\widetilde{\mathcal{V}}_{2m})=0,$
(53)

where $\mathcal{R}_{m}$ is the tensor residual given by
$\mathcal{R}_{m}=\mathcal{X}_{m}-{\cal A}\ast\mathcal{X}_{m}\ast{\cal
A}^{T}-{\cal B}\ast{\cal B}^{T}$. By developing the condition described in
(53), we get the low order discrete Lyapunov matrix equation verified by
$Y_{2m}$, which can be expressed as follows

$Y_{2m}-T_{m}Y_{2m}T_{m}^{T}=B_{m}B_{m}^{T},$ (54)

where $T_{m}\in\mathbb{R}^{2m\times 2m}$ is a Hessenberg matrix obtained from
Algorithm 2 applied to the pair (${\cal A},{\cal B}$) and
$B_{m}=\widetilde{\mathcal{V}}_{2m}^{T}\diamond{\cal B}$. A simplification to
$B_{m}$ goes as

$B_{m}=\widetilde{\mathcal{V}}_{2m}^{T}\diamond{\cal
B}=\omega_{11}e_{1}^{(2m)},\quad\text{since}\,\,{\cal
B}=\omega_{11}\mathcal{V}_{1}^{1}\,\,(\text{from Algorithm
\ref{tenextgloalgo}})$

where $\mathcal{V}_{1}^{1}$ is defined using the MATLAB notation as
$\mathcal{V}_{1}^{1}=\mathcal{V}_{1}(:,:,:,1:K_{2})$, $\omega_{11}$ is a
scalar defined in Step $2$ of Algorithm 2 and $e_{1}^{(2m)}$ is the first
column of the identity matrix $I_{2m}\in\mathbb{R}^{2m\times 2m}$. The
following result shows how to compute the residual in an efficient way.

######  Theorem 1.

Assume that $m$ iterations of Algorithm 2 have been run. Then

$\|\mathcal{R}_{m}\|\leq r_{m},$ (55)

where

$r_{m}=\sqrt{2\|T_{m}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\|_{F}^{2}+\|(Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T})^{T}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\|_{F}^{2}}.$

$E_{(2m)}^{m}$ is the last $2m\times 2$ column block of the identity matrix
$I_{2m}$.

###### Proof.

By using the definition of the approximation $\mathcal{X}_{m}$ and the
decompositions that we get from Algorithm 2, we end up with the following

$\begin{split}\mathcal{R}_{m}&=\mathcal{X}_{m}-{\cal
A}\ast\mathcal{X}_{m}\ast{\cal A}^{T}-{\cal B}\ast{\cal B}^{T}\\\
&=\widetilde{\mathcal{V}}_{2m}\ast(Y_{2m}\circledast{\cal
I}_{2m})\ast\widetilde{\mathcal{V}}_{2m}^{T}-\underbrace{{\cal
A}\ast(\widetilde{\mathcal{V}}_{2m}\ast(Y_{2m}\circledast{\cal
I}_{2m})\ast\widetilde{\mathcal{V}}_{2m}^{T})\ast{\cal
A}^{T}}_{(1)}-\underbrace{{\cal B}\ast{\cal B}^{T}}_{(2)}.\end{split}$ (56)

The term $(1)$ could be simplified as follows

$\begin{array}[]{lll}{\cal
A}\ast(\widetilde{\mathcal{V}}_{2m}\ast(Y_{2m}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T})\ast{\cal
A}^{T}&=&(\widetilde{\mathcal{V}}_{2m}\ast(T_{m}Y_{2m}\circledast{\cal
I}_{K})+\mathcal{V}_{2(m+1)}\ast(T_{m+1,m}(E_{(2m)}^{m})^{T}Y_{2m}\circledast{\cal
I}_{K}))\\\ &\ast&((T_{m}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}+(E_{(2m)}^{m}T_{m+1,m}^{T}\circledast{\cal
I}_{K})\ast\mathcal{V}_{2(m+1)}^{T})\\\
&=&\widetilde{\mathcal{V}}_{2m}\ast(T_{m}Y_{2m}T_{m}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}\\\
&+&\widetilde{\mathcal{V}}_{2m}\ast(T_{m}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\circledast{\cal
I}_{K})\ast\mathcal{V}_{2(m+1)}^{T}\\\
&+&\mathcal{V}_{2(m+1)}\ast(T_{m+1,m}(E_{(2m)}^{m})^{T}Y_{2m}T_{m}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}\\\
&+&\mathcal{V}_{2(m+1)}\ast(T_{m+1,m}(E_{(2m)}^{m})^{T}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\circledast{\cal
I}_{K})\ast\mathcal{V}_{2(m+1)}^{T}.\end{array}$

From Algorithm 2, we know that ${\cal B}=\omega_{1,1}\mathcal{V}_{1}^{1}$,
where $\mathcal{V}_{1}^{1}$ is defined using the MATLAB notation as
$\mathcal{V}_{1}^{1}=\mathcal{V}_{1}(:,:,:,1:K_{2})$ and $\omega_{1,1}$ is a
scalar defined in Step $2$ of Algorithm 2, then the term $(2)$ could be
written as

${\cal B}\ast{\cal
B}^{T}=(\omega_{1,1})^{2}\widetilde{\mathcal{V}}_{2m}\ast(e_{1}^{(2m)}(e_{1}^{(2m)})^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}.$

Gathering all these results, and by considering
$Y_{2m}-T_{m}Y_{2m}T_{m}^{T}-\omega_{1,1}^{2}(e_{1}^{(2m)}(e_{1}^{(2m)})^{T})=0$
mentioned in (54), we end up with

$\mathcal{R}_{m}=\widetilde{\mathcal{V}}_{2(m+1)}\ast\left[\begin{pmatrix}0&-T_{m}Y_{2m}E_{(2m)}^{m}T_{m+1,m}\\\
-T_{m+1,m}^{T}(E_{(2m)}^{m})^{T}Y_{2m}T_{m}^{T}&T_{m+1,m}(E_{(2m)}^{m})^{T}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\end{pmatrix}\circledast{\cal
I}_{K}\right]\ast\widetilde{\mathcal{V}}_{2(m+1)}^{T}.$

It follows that

$\|\mathcal{R}_{m}\|\leq\sqrt{2\|T_{m}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\|_{F}^{2}+\|(Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T})^{T}Y_{2m}E_{(2m)}^{m}T_{m+1,m}^{T}\|_{F}^{2}},$

which ends the proof. ∎

######  Remark 3.

The approximation $\mathcal{X}_{m}$ could be expressed in a factored form.
This allows us to compute $\mathcal{X}_{m}$ in a efficient way. We start by
employing a Singular Value Decomposition (SVD) to $Y_{2m}$, i.e.,
$Y_{2m}=P\Sigma Q$, then we consider some tolerance dtol and define $P_{r}$,
$Q_{r}$ as the first $r$ columns of $P$ and $Q$ corresponding respectively to
the $r$ singular values of magnitude greater than dtol. By setting
$\Sigma_{r}=\text{diag}(\sigma_{1}\ldots,\sigma_{r})$, we get $Y_{2m}\approx
P_{r}\Sigma_{r}Q_{r}^{T}$ and the factorization goes as

$\mathcal{X}_{m}\approx\widetilde{\mathcal{V}}_{2m}\ast(P_{r}\Sigma_{r}Q_{r}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}\approx{\cal Z}_{1}\ast{\cal
Z}_{2}^{T},$ (57)

where ${\cal
Z}_{1}=\widetilde{\mathcal{V}}_{2m}\ast(P_{r}(\Sigma_{r})^{1/2}\circledast{\cal
I}_{K})$ and ${\cal Z}_{2}=((\Sigma_{r})^{1/2}Q_{r}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}$. The same process could be
followed to get the approximation in a factored form using the global extended
tensor Krylov subspace, in this case ${\cal
Z}_{1}=\widetilde{\mathcal{V}}_{2m}\ast(P_{r}(\Sigma_{r})^{1/2}\circledast{\cal
I}_{K})$ and ${\cal Z}_{2}=((\Sigma_{r})^{1/2}Q_{r}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}$.

All the results described in this section enable us to provide with the
following algorithm

Algorithm 3 The Lyapunov tensor extended global Arnoldi algorithm (TLTEGAA)

  1. 1.

Inputs: ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$,
${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$,
tolerances $\epsilon$, dtol and a number $m_{max}$ of maximum iterations.

  2. 2.

Outputs: the approximate solution $\mathcal{X}_{m}\approx{\cal Z}_{1}\ast{\cal
Z}_{2}^{T}$.

  3. 3.

For $m=1,\cdots,m_{max}$

  4. 4.

Use Algorithm 2 to built $\widetilde{\mathcal{V}}_{2m}$ an orthonormal tensor
associated to the tensor extended Krylov subspace (25) and
$H_{m}\in\mathbb{R}^{2m\times 2m}$ a Hessenberg matrix. We use Proposition 8
to compute the Hessenberg matrix $T_{m}\in\mathbb{R}^{2m\times 2m}$.

  5. 5.

Solve the low-dimensional disctete Lyapunov equation (54) using the MATLAB
function dlyap.

  6. 6.

Compute the residual norm $\|\mathcal{R}_{m}\|$ using Theorem 1, and if it is
less than $\epsilon$, then

    1. (a)

compute the SVD of $Y_{2m}=P\Sigma Q^{T}$ where
$\Sigma=diag[\sigma_{1},\cdots,\sigma_{2m}],$

    2. (b)

determine $r$ such that $\sigma_{r+1}<{\tt dtol}\leq\sigma_{r}$, set
$\Sigma_{r}=diag[\sigma_{1},\cdots,\sigma_{r}]$ and compute ${\cal
Z}_{1}=\widetilde{\mathcal{V}}_{2m}\ast(P_{r}(\Sigma_{r})^{1/2}\circledast{\cal
I}_{K})$ and ${\cal Z}_{2}=((\Sigma_{r})^{1/2}Q_{r}^{T}\circledast{\cal
I}_{K})\ast\widetilde{\mathcal{V}}_{2m}^{T}$,

end if.

  7. 7.

End For

###  5.2 Tensor extended block Arnoldi process for Lyapunov equations

In this section, we define a tensor extended block Krylov subspace via
Einstein product that will be used to solve the tensor discrete Lyapunov
equation (50). A description of the classic block Krylov subspace methods via
Einstein product has been performed in [17]. For ${\cal
A}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times J_{1}\times\ldots\times
J_{N}}$ and ${\cal B}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}}$, the $2m$-th block extended tensor Krylov
subspace associated to the pair $({\cal A},{\cal B})$ is defined by

${\cal K}_{m}^{b-ext}({\cal A},{\cal B})=\text{Range}\left\\{{\cal A}^{-\star
m}\ast{\cal B},\ldots,{\cal A}^{-1}\ast{\cal B},{\cal B},\ldots,{\cal
A}^{\star\,m-1}\ast{\cal B}\right\\},$ (58)

where ${\cal A}^{\star k}=\underbrace{{\cal A}\ast\ldots\ast{\cal A}}_{k}$,
${\cal A}^{\star-k}=\underbrace{{\cal A}^{-1}\ast\ldots\ast{\cal A}^{-1}}_{k}$
and ${\cal A}^{0}={\cal I}$ is the identity tensor. Next, we define the notion
of upper/lower triangular tensors which will be used later on.

######  Definition 11.

Let ${\cal U}$ and ${\cal L}$ be two tensors in the space
$\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times I_{1}\times\ldots\times
I_{N}}$

  * -

${\cal U}$ si an upper triangular tensor if the entries
$u_{j_{1},\ldots,j_{N},i_{1},\ldots,i_{N}}=0$ when $i\text{vec}(j,J)\geq
i\text{vec}(j,J)$.

  * -

${\cal L}$ si a lower triangular tensor if the entries
$r_{j_{1},\ldots,j_{N},i_{1},\ldots,i_{N}}=0$ when $i\text{vec}(j,J)\leq
i\text{vec}(j,J)$,

where $i\text{vec}(\cdot,\cdot)$ is the mapping index mentioned in Definition
4.

Analogously to the QR decomposition in the matrix case [15], a similar
definition to such decomposition of a tensor $\mathcal{P}$ in the space
$\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times K_{1}\times\ldots\times
K_{M}}$ is defined as follows

$\mathcal{P}={\cal Q}\ast\mathcal{R},$ (59)

where ${\cal Q}\in\mathbb{R}^{J_{1}\times\ldots\times J_{N}\times
K_{1}\times\ldots\times K_{M}}$ is orthogonal i.e., ${\cal Q}^{T}\ast{\cal
Q}={\cal I}_{K}$ and $\mathcal{R}\in\mathbb{R}^{K_{1}\times\ldots\times
K_{M}\times K_{1}\times\ldots\times K_{M}}$ is an upper triangular tensor
(Definition 11). Similarly to the matrix extended block Krylov Arnoldi
algorithm [38], an analogous algorithm in a tensor format, particularly for
the $4th$-order tensor which is the main subject of this paper as mentioned
previously, could be generated as follows

Algorithm 4 Tensor extended block Arnoldi algorithm

  1. 1.

Input : ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$,
${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$ and a
fixed integer $m$.

  2. 2.

Output : $\widetilde{\mathcal{V}}_{2(m+1)}\in\mathbb{R}^{J_{1}\times
J_{2}\times K_{1}\times 2(m+1)K_{2}}$ and
$\overline{H}_{m}\in\mathbb{R}^{K_{1}\times 2(m+1)K_{2}\times K_{1}\times
2mK_{2}}$.

  3. 3.

Compute $QR$-decomp. to $|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times 2K_{2}}$ i.e.,
$|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2}=\mathcal{V}_{1}\ast\mathcal{R}_{0}$. (MATLAB notation
$\widetilde{\mathcal{V}}_{2(m+1)}(:,:,:,1:2K_{2})=\mathcal{V}_{1}$)

  4. 4.

for $j=1\cdots m$

  5. 5.

Set $\mathcal{W}=|{\cal A}\ast\mathcal{V}_{j}^{1},{\cal
A}^{-1}\ast\mathcal{V}_{j}^{2}|_{2}$, where
$\left\\{\begin{array}[]{lll}\mathcal{V}_{j}^{1}&=&\mathcal{V}_{j}(:,:,:,(j-1)K_{2}+1:jK_{2}),\\\
\mathcal{V}_{j}^{2}&=&\mathcal{V}_{j}(:,:,:,jK_{2}+1:(j+1)K_{2}),\end{array}\right.$

  6. 6.

for $i=1,\ldots,j$.

  7. 7.

${\cal H}_{i,j}=\mathcal{V}_{i}\ast\mathcal{W}$.

  8. 8.

$\mathcal{W}=\mathcal{W}-\mathcal{V}_{i}\ast{\cal H}_{i,j}$.

  9. 9.

end for.

  10. 10.

Compute $QR$-decomp. to $\mathcal{W}$, i.e.,
$\mathcal{W}=[\mathcal{V}_{j+1},{\cal H}_{j+1,j}].$

  11. 11.

end for.

We start our process by applying Algorithm 4 defined above to the pair $({\cal
A},{\cal B})$. After $m$ iterations, the following decomposition could be
obtained

$\begin{array}[]{lll}{\cal
A}\ast\widetilde{\mathcal{V}}_{2m}&=\widetilde{\mathcal{V}}_{2(m+1)}\ast\overline{{\cal
T}}_{m}\\\ &=\widetilde{\mathcal{V}}_{2m}\ast{\cal
T}_{m}+\mathcal{V}_{2(m+1)}\ast({\cal T}_{m+1,m}\ast{\cal
E}_{2m}^{T}),\end{array}$ (60)

where $\widetilde{\mathcal{V}}_{2m}$ is the tensor basis associated with
${\cal K}_{m}^{b-ext}({\cal A},{\cal B})$. Tensors $\overline{{\cal T}}_{m}$
and ${\cal T}_{m}$ are defined as follows

$\overline{{\cal T}}_{m}=\widetilde{\mathcal{V}}_{2(m+1)}^{T}\ast({\cal
A}\ast\widetilde{\mathcal{V}}_{2m})\in\mathbb{R}^{K_{1}\times
2(m+1)K_{2}\times K_{1}\times 2mK_{2}},\quad{\cal
T}_{m}=\widetilde{\mathcal{V}}_{2m}^{T}\ast({\cal
A}\ast\widetilde{\mathcal{V}}_{2m})\in\mathbb{R}^{K_{1}\times 2mK_{2}\times
K_{1}\times 2mK_{2}}.$

Finally the tensor ${\cal E}_{2m}$ is get it from the identity tensor ${\cal
I}_{2m}\in\mathbb{R}^{K_{1}\times 2mK_{2}\times K_{1}\times 2mK_{2}}$ as
${\cal E}_{2m}={\cal I}_{2m}(:,:,:,2(m-1)K_{2}+1:2mK_{2})$. Analogously to the
matrix case, the tensor ${\cal T}_{m}$ can be defined as the restriction of
the tensor ${\cal A}$ in the extended block tensor Krylov subspace ${\cal
K}_{m}^{b-ext}({\cal A},{\cal B})$. Notice that the tensor $\overline{{\cal
T}}_{m}$ can be expressed as

$\overline{{\cal T}}_{m}=\begin{vmatrix}{\cal T}_{1,1}&{\cal
T}_{1,2}&\ldots&{\cal T}_{1,m}\\\ {\cal T}_{2,1}&{\cal T}_{2,2}&\ldots&{\cal
T}_{2,m}\\\ \vdots&\ddots&\ddots&\vdots\\\ \mathcal{O}&\ldots&{\cal
T}_{m,m-1}&{\cal T}_{m,m}\\\ \mathcal{O}&\ldots&\mathcal{O}&{\cal
T}_{m+1,m}\end{vmatrix}\in\mathbb{R}^{K_{1}\times 2(m+1)K_{2}\times
K_{1}\times 2mK_{2}},$ (61)

where the notation $|\cdot|$ is the definition of block tensors given in
Section 1.2. The tensor ${\cal T}_{m}$ is obtained from $\overline{{\cal
T}}_{m}$ by deleting the last row block tensor
$\big{|}\mathcal{O}\quad\mathcal{O}\quad\ldots\quad{\cal H}_{m+1,m}\big{|}$.

######  Remark 4.

Generally, the tensor $\overline{{\cal H}}_{m}\in\mathbb{R}^{K_{1}\times
2(m+1)K_{2}\times K_{1}\times 2mK_{2}}$ is of small dimension, then the
computation of the matricization of such tensor denoted by
$mat(\overline{{\cal H}}_{m})=\overline{H}_{m}$ is feasible. We denote by
$\overline{T}_{m}$ the matrix that we get from the matricization of tensor
$\overline{{\cal T}}_{m}$. We recall that $\overline{T}_{m}$ could be computed
efficiently and directly using $\overline{H}_{m}$ as mentioned in [37] without
requiring any extra product. The purpose behind this process is to build
${\cal T}_{m}$ needed afterwards without requiring any extra tensor product
with the large scale tensor ${\cal A}$.

We seek for an approximate solution $\mathcal{X}_{m}$ to $\mathcal{X}$
fulfilling the discrete Lyapunov equation (51). Consider the following form of
the approximation

$\mathcal{X}_{m}=\widetilde{\mathcal{V}}_{2m}\ast\mathcal{Y}_{2m}\ast\widetilde{\mathcal{V}}_{2m}^{T},$
(62)

where $\widetilde{\mathcal{V}}_{2m}$ is the tensor basis associated with the
tensor Krylov subspace obtained by using Algorithm 4 to the pair $({\cal
A},{\cal B})$. We obtain the tensor
$\mathcal{Y}_{2m}\in\mathbb{R}^{K_{1}\times 2mK_{2}\times K_{1}\times
2mK_{2}}$ via the following condition

$\widetilde{\mathcal{V}}_{2m}^{T}\ast(\mathcal{R}_{m}\ast\widetilde{\mathcal{V}}_{2m})=0,$
(63)

and $\mathcal{R}_{m}$ is the residual tensor given by
$\mathcal{R}_{m}=\mathcal{X}_{m}-{\cal A}\ast\mathcal{X}_{m}\ast{\cal
A}^{T}-{\cal B}\ast{\cal B}^{T}$. We develop the condition (63) using the
decomposition (60) together with orthogonality condition (i.e.,
$\widetilde{\mathcal{V}}_{2m}^{T}\ast\widetilde{\mathcal{V}}_{2m}={\cal
I}_{2m}$), we find that the tensor $\mathcal{Y}_{2m}$ fulfil the following low
dimensional Lyapunov tensor equation

$\mathcal{Y}_{m}-{\cal T}_{m}\ast\mathcal{Y}_{m}\ast{\cal T}_{m}^{T}={\cal
B}_{m}\ast{\cal B}_{m}^{T}.$ (64)

The right hand side could be simplified to the following ${\cal
B}_{m}\ast{\cal B}_{m}^{T}$. We know from Algorithm 4 that

$\begin{split}|{\cal B},\,{\cal A}^{-1}\ast{\cal
B}|_{2}&=\mathcal{V}_{1}\ast\mathcal{R}_{0}\\\
&=|\mathcal{V}_{1}^{1},\mathcal{V}_{1}^{2}|_{2}\ast\begin{vmatrix}\mathcal{R}_{0}^{(1,1)}&\mathcal{R}_{0}^{(1,2)}\\\
\mathcal{O}&\mathcal{R}_{0}^{(2,2)}\end{vmatrix},\end{split}$

where
$\mathcal{V}_{1}^{i}=\mathcal{V}_{1}(:,:,:,(i-1)K_{2}+1:iK_{2})\in\mathbb{R}^{K_{1}\times
K_{2}\times K_{1}\times K_{2}}$,
$\mathcal{R}_{0}^{(i,j)}=\mathcal{R}_{0}(:,(i-1)K_{2}+1:iK_{2},:,(j-1)K_{2}+1:jK_{2})\in\mathbb{R}^{K_{1}\times
K_{2}\times K_{1}\times K_{2}}$. The notation $|\cdot|_{2}$ is the definition
of row block tensor given in Definition 6 and the upper triangular tensor
$\mathcal{R}_{0}$ defined in Definition 11 can be presented as described above
using notation from Section 1.2. Using Proposition 3, we can conclude the
following

${\cal B}=\mathcal{V}_{1}^{1}\ast\mathcal{R}_{0}^{(1,1)},$

thus

$\begin{split}{\cal B}_{m}&=\widetilde{\mathcal{V}}_{2m}^{T}\ast{\cal B}\\\
&={\cal E}_{2m}^{1}\ast\mathcal{R}_{0}^{(1,1)},\end{split}$ (65)

where ${\cal E}_{2m}^{1}$ can be extracted from the identity tensor ${\cal
I}_{2m}\in\mathbb{R}^{K_{1}\times 2mK_{2}\times K_{1}\times 2mK_{2}}$ as
follows ${\cal E}_{2m}^{1}={\cal I}_{2m}(:,:,:,1:K_{2}).$ In the next result,
we show how to simplify the expression of the residual $\mathcal{R}_{m}$, this
yields to a proper computation of the residual tensorial norm
$\|\mathcal{R}_{m}\|$. Our process is stopped whenever the residual norm
$\|\mathcal{R}_{m}\|$ is lower than a chosen tolerance $\epsilon$
(i.e.$\|\mathcal{R}_{m}\|<\epsilon$ ).

######  Theorem 2.

Assume that $m$ iterations have been run using Algorithm 4, then

$\|\mathcal{R}_{m}\|=r_{m},$ (66)

where

$r_{m}=\sqrt{2*\|{\cal T}_{m}\ast\mathcal{Y}_{2m}\ast{\cal E}_{2m}\ast{\cal
T}_{m+1,m}^{T}\|^{2}+\|{\cal T}_{m+1,m}\ast{\cal
E}_{2m}^{T}\ast\mathcal{Y}_{2m}\ast{\cal E}_{2m}\ast{\cal
T}_{m+1,m}^{T}\|^{2}},$

and $\|\cdot\|$ is the Frobenius tensor norm defined in Definition 3.

###### Proof.

We plug $\mathcal{X}_{m}$ in the expression of $\mathcal{R}_{m}$ and by using
decomposition (60) and the fact that $\mathcal{Y}_{2m}-{\cal
T}_{m}\ast\mathcal{Y}_{2m}\ast{\cal T}_{m}^{T}-{\cal
E}_{2m}^{1}\ast\mathcal{R}_{0}^{(1,1)}\ast({\cal
E}_{2m}^{1}\ast\mathcal{R}_{0}^{(1,1)})^{T}=0$ (i.e., (64)), we get

$\mathcal{R}_{m}=\widetilde{\mathcal{V}}_{2(m+1)}\ast\left[\begin{vmatrix}0&-{\cal
T}_{m}\ast\mathcal{Y}_{2m}\ast{\cal E}_{2m}\ast{\cal T}_{m+1,m}^{T}\\\ -({\cal
T}_{m}\ast\mathcal{Y}_{m}\ast{\cal E}_{2m}\ast{\cal T}_{m+1,m}^{T})^{T}&{\cal
T}_{m+1,m}\ast{\cal E}_{2m}^{T}\ast\mathcal{Y}_{2m}\ast{\cal E}_{2m}\ast{\cal
T}_{m+1,m}^{T}\end{vmatrix}\right]\ast\widetilde{\mathcal{V}}_{2(m+1)}^{T}.$

Then we obtain the desired result

$\|\mathcal{R}_{m}\|=\sqrt{2*\|{\cal T}_{m}\ast\mathcal{Y}_{2m}\ast{\cal
E}_{2m}\ast{\cal T}_{m+1,m}^{T}\|^{2}+\|{\cal T}_{m+1,m}\ast{\cal
E}_{2m}^{T}\ast\mathcal{Y}_{2m}\ast{\cal E}_{2m}\ast{\cal
T}_{m+1,m}^{T}\|^{2}}.$

∎

Thanks to all the results described above, we are now able to present the
following algorithm

Algorithm 5 The Lyapunov tensor extended block Arnoldi algorithm (TLTEBAA)

  1. 1.

Inputs: ${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$,
${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$,
tolerances $\epsilon$, dtol and a number $m_{max}$ of maximum iterations.

  2. 2.

Outputs: the approximate solution $\mathcal{X}_{m}\approx{\cal Z}_{1}\ast{\cal
Z}_{2}^{T}$.

  3. 3.

For $m=1,\cdots,m_{max}$

  4. 4.

Use Algorithm 4 to built $\widetilde{\mathcal{V}}_{2m}$ an orthonormal tensor
associated to the tensor Krylov subspace ${\cal K}_{m}^{b-ext}({\cal A},{\cal
B})$ and ${\cal H}_{m}\in\mathbb{R}^{K_{1}\times 2mK_{2}\times K_{1}\times
2mK_{2}}$ a type Hessenberg form tensor.

  5. 5.

Solving the low-dimensional Lyapunov tensor equation (64) after matricization
using the MATLAB function dlyap.

  6. 6.

Compute the residual norm $\|\mathcal{R}_{m}^{e}\|$ using Theorem 2, and if it
is less than $\epsilon$, then

    1. (a)

Compute the SVD of $\mathcal{Y}_{m}=\mathcal{P}\ast\Sigma\ast{\cal Q}^{T}$
where $\Sigma$ is a diagonal tensor contains singular values
$\left\\{\sigma_{1},\cdots,\sigma_{2m}\right\\}.$ This definition of tensor
SVD via Einstein product could be found in [9].

    2. (b)

Determine $r$ such that $\sigma_{r+1}<{\tt dtol}\leq\sigma_{r}$, set
$\Sigma_{r}$ as a diagonal tensor with only the first $r$ singular values,
then we compute ${\cal
Z}_{1}=\widetilde{\mathcal{V}}_{2m}\ast(\mathcal{P}_{r}\ast(\Sigma_{r})^{1/2})$
and ${\cal Z}_{2}=((\Sigma_{r})^{1/2}\ast{\cal
Q}_{r}^{T})\ast\widetilde{\mathcal{V}}_{2m}^{T}$. $\mathcal{P}_{r}$ and ${\cal
Q}_{r}$ are defined as $\mathcal{P}_{r}=\mathcal{P}(:,:,:,1:2rK_{2})$ and
${\cal Q}_{r}={\cal Q}(:,:,:,1:2rK_{2}).$

end if.

  7. 7.

End For

######  Remark 5.

Dealing with the continuous time MLTI systems (11) lead to the appearance of
continuous time Lyapunov equations of the form ${\cal
A}\ast\mathcal{X}-\mathcal{X}\ast{\cal A}^{T}={\cal B}\ast{\cal B}^{T}$. The
methods described above for finding an efficient approximation to the
discrete-time Lyapunov equations (50) can also be applied to the continuous-
time Lyapunov equations.

##  6 Numerical experiments

In this section, numerical examples are performed to validate the efficiency
of the algorithms described above based on the tensor extended global and
block Krylov subspaces. The numerical results were obtained using MATLAB
R2018a on a computer with Intel ${}^{\text{\textregistered}}$ core i7 at
2.6GHz and 16Gb of RAM. We need to mention that all the algorithms described
here, have been implemented based on the MATLAB tensor toolbox established by
Kolda et al., [26].
We start our numerical tests by showing the results obtained by the tensor
extended global algorithm to get a reduced MLTI system from the original one.
The tensorial structure of the constructed reduced MLTI system is described in
(41). We are interested in the input-output behaviour of the reduced MLTI
system. Such behaviour could be measured by the closeness of the original and
reduced transfer functions, denoted by ${\cal F}(s)$ and ${\cal F}_{m}(s)$,
respectively. Analogously to the matrix case ([10]), we can define the
following norm named $h_{\infty}$-norm (cf. [2]) to describe the magnitudes
needed for plotting the two transfer functions and see how close the reduced
MLTI system to the original one

$\|{\cal F}(.)\|_{\infty}=\|\Psi({\cal
F}(.))\|_{\infty}=\sup_{\theta\in[0,2\pi]}\sigma_{max}(\Psi({\cal
F}(e^{i\theta}))).$ (67)

Example 1. We consider the following data

  * -

${\cal A}\in\mathbb{R}^{N\times N\times N\times N}$, with $N=100$. Here,
${\cal A}$ is constructed from a tensorization of a triangular matrix
$A\in\mathbb{R}^{N\,N\times N\,N}$ constructed using spdiags MATLAB function.

  * -

${\cal B},{\cal C}^{T}\in\mathbb{R}^{N\times N\times K_{1}\times K_{2}}$ are
chosen as a sparse and random tensors respectively with $K_{1}=3$ and
$K_{2}=5$.

The choice of $m=5$ which is the dimension of the tensor extended global
Krylov subspace where we project the original MLTI system was sufficient to
obtain the results described in Figure 1.

Fig. 1: The frequency responses of the original and reduced MLTI systems
(left) and the error norms $\|{\cal F}-{\cal F}_{m}\|_{\infty}$ (right).

Example 2. This example is from ([12]) and it concerns the evolution of heat
distribution in a solid medium over time. The partial differential equation
that describes this evolution is given by the following equation and named the
$2D$ heat equation

$\left\\{\begin{split}\dfrac{\partial}{\partial
t}\phi(t,x)&=c^{2}\dfrac{\partial^{2}}{\partial
x^{2}}\phi(t,x)+\delta(x)u_{t},\quad x\in D\,(\text{square
}D=[-\pi^{2},\pi]),\\\ \phi(t,x)&=0,\quad x\in\partial D,\end{split}\right.$
(68)

where $c>0$, $u_{t}$ is a one dimensional control input and $\delta(x)$ is the
Dirac delta function centred at zeros. According to [12], a second-order
central difference is used to approximate the Laplacien and a first-order
difference in time, this leading to a MLTI system of the form (10),where
$\mathcal{X}_{k}\in\mathbb{R}^{N\times N}$ called the 2D-temperature field at
instance $k$. For further information, we refer the readers to [12], the same
problem has been treated in a different way in [9]. We consider the following
system tensors

  * -

The tensor ${\cal A}\in\mathbb{R}^{N\times N\times N\times N}$ with $N=128$ is
the the tensorization of $\frac{c^{2}\Delta
t}{h^{2}}\Delta_{dd}\in\mathbb{R}^{N^{2}\times N^{2}}$, where $\Delta_{dd}$ is
the discrete Laplacian on a rectangular grid with a Dirichlet boundary
condition (see [12] for more details).

  * -

The constructed tensor ${\cal B}$ in [12] is of dimension $N\times N\times
1\times 1$. Here, we change ${\cal B}$ and we choose it as a sparse tensor in
$\mathbb{R}^{N\times N\times K_{1}\times K_{2}}$ with $K_{1}=3$ and $K_{2}=5$.

  * -

To complete the system, we will need the output tensor ${\cal C}$. We choose
it as random tensor in $\mathbb{R}^{K_{1}\times K_{2}\times N\times N}$ with
$K_{1}=3$ and $K_{2}=5$.

  * -

We set $m=20$, the associated dimension of the extended tensor Krylov subspace
(25).

Fig. 2: The frequency responses of the original and reduced MLTI systems
(left) and the error norms $\|{\cal F}-{\cal F}_{m}\|_{\infty}$ (right).

Example 3. For this experiment, we show the results obtained from the tensor
Balanced Truncation method explained in Section 4, we consider the following
data

  * -

${\cal A}\in\mathbb{R}^{N\times N\times N\times N}$, with $N=100$. Here,
${\cal A}$ is constructed from a tensorization of a matrix
$A\in\mathbb{R}^{N^{2}\times N^{2}}$ constructed using the identity matrix
speye perturbed by sprandn MATLAB function.

  * -

${\cal B},{\cal C}^{T}\in\mathbb{R}^{N\times N\times K_{1}\times K_{2}}$ are
chosen as a random tensors with $K_{1}=5$ and $K_{2}=5$.

  * -

We use Algorithm 3 to solve the two discrete Lyapunov equations (45) by
setting $m=20$ (i.e., max number of iterations) and the tolerance
$\epsilon=1e^{-6}$.

Fig. 3: The frequency responses of the original and reduced MLTI systems
(left) and the error norms $\|{\cal F}-{\cal F}_{m}\|_{\infty}$ (right).

For the next experiment, we show the results obtained by the methods described
in Algorithms 4 and 5. These methods are seeking for an efficient approximate
solution of the discrete Lyapunov tensor equation $\mathcal{X}-{\cal
A}\ast\mathcal{X}\ast{\cal A}^{T}={\cal B}\ast{\cal B}^{T}$. We stopped our
iterations whenever $\|\mathcal{R}_{m}\|<\epsilon$, where $\|\cdot\|$ is the
Frobenius norm associated to tensors defined in Section 1 and $\epsilon$ is
some tolerance that will be mentioned in each example. We compared our methods
based on the tensor extended global and block Krylov subspaces, in Algorithm 3
(for the global process) and in Algorithm 5 (for the block process) with the
classic global Arnoldi described in Algorithm 1 and block Arnoldi process
described in [17]. Consider the following data

  * -

${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$, with
$J_{1}=J_{2}$ in $\left\\{32,64,128\right\\}.$ Here, ${\cal A}$ is a random
sparse tensor (i.e., MATLAB func. sptenrand).

  * -

${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$ is chosen
as a sparse tensor with $K_{1}=3$ and $K_{2}=5$.

  * -

Th tolerance $\epsilon$ is set to $10^{-6}$.

Fig. 4: Consumed time by the four methods to reach the criterion convergence
(left) and the residual norm (right).

In the left plot of Figure 4, we show the cpu-time consumed by the four
methods to achieve the convergence criterion (i.e.,
$\|\mathcal{R}_{m}\|<1e^{-6}$). The dimension of the tensor extended Krylov
subspace was fixed to $m=20$, and since this latter is of $2m$ dimension, we
used $m_{1}=2m$ the dimension of the tensor classic Krylov subspace in order
to realise a fair comparison. This latter is considered for all the examples
described here. The tensor extended block and global Krylov subspace methods
have successfully achieved the desired result in a few seconds as mentioned in
the right plot of Figure 4, contrary to the tensor classic block and global
Krylov subspace methods where no convergence has been guaranteed and more
iterations were need.

For the next experiment, the tolerance was set to $\epsilon=10^{-8}$ and the
following data has been considered

  * -

${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$, we set
$J_{1}=J_{2}$ in $\left\\{32,64,128,256\right\\}.$ Here, ${\cal A}$ is
constructed from a tensorization of a triangular matrix
$A\in\mathbb{R}^{J_{1}J_{2}\times J_{1}J_{2}}$ constructed using spdiags
MATLAB function.

  * -

${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$ is chosen
as a random tensor with $K_{1}=3$ and $K_{2}=5$. In Figure, ${\cal B}$ is a
random tensor with $K_{1}=5$ and $K_{2}=6$.

Fig. 5: Consumed time by the four methods to reach the criterion convergence.

In Figure 5, we represent the computation time used for the four methods to
reach the convergence criterion. In Figure 6, we plotted the cpu-time needed
to achieve a good convergence using two examples of ${\cal B}$; for the left
plot we chose ${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times 5\times 6}$, and
for the right plot, we show the results depicted from choosing ${\cal
B}\in\mathbb{R}^{J_{1}\times J_{2}\times 6\times 8}$.

Fig. 6: Consumed time by the four methods to reach the criterion convergence,
where the right hand side ${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times
5\times 6}$ (left) ${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times 6\times 8}$
(right).

In Table 1, we present the results obtained from the tensor classic and
extended Krylov subspace methods using the block and global process to solve
high order discrete Lyapunov tensor equations. We stopped our iterations when
the residual norm is less than $\epsilon=10^{-8}$. We considered the following
data

  * -

${\cal A}\in\mathbb{R}^{J_{1}\times J_{2}\times J_{1}\times J_{2}}$, we set
$J_{1}=J_{2}$ in $\left\\{512,1024\right\\}.$ Here, ${\cal A}$ is constructed
from a tensorization of a triangular matrix $A\in\mathbb{R}^{J_{1}J_{2}\times
J_{1}J_{2}}$ constructed using spdiags MATLAB function. After that we use
spreshape to get the tensor form.

  * -

${\cal B}\in\mathbb{R}^{J_{1}\times J_{2}\times K_{1}\times K_{2}}$ is chosen
as a random tensor with $K_{1}=5$ and $K_{2}=6$.

Table 1: Results gained from the four methods for solving high order Lyapunov
tensor equations.

Size of ${\cal A}$ ($J_{1}\times J_{2}\times J_{1}\times J_{2}$) | Classic Kry. Bloc | Classic Kry. Global | Extended Kry. Bloc | Extended Kry. Global
---|---|---|---|---
| cpu-time Resi-norm Iter. | cpu-time Resi-norm Iter. | cpu-time Resi-norm Iter. | cpu-time Resi-norm Iter.
$J_{1}=J_{2}=512$ |  42.84 4.07 $\times$ 10-09 25 |  36.25 5.41$\times$ 10-09 26 |  9.63 2.81$\times$ 10-09 7 |  14.37 1.83$\times$ 10-10 8
$J_{1}=J_{2}=1024$ |  172.02 9.42 $\times$ 10-09 26 |  149.88 6.78$\times$ 10-09 27 |  48.91 1.75$\times$ 10-10 8 |  57.14 7.10$\times$ 10-10 8

## Conclusion

In this paper, we explored the application of multilinear algebra in reducing
the order of multidimentional linear time invariant (MLTI) systems. We used
tensor Krylov subspace methods as key tools, which involve approximating the
system solution within a low-dimensional subspace. We introduced the tensor
extended block and global Krylov subspaces and the corresponding Arnoldi based
processes. Using these methods, we developed a model reduction using
projection techniques. We have also shown how these methods could be used to
solve large-scale Lyapunov tensor equations that are required in the balanced
truncation method which is a technique for order reduction. We demonstrated
how to extract approximate solutions via the Einstein product using the tensor
extended block Arnoldi and the extended global Arnoldi processes.

## References

  * [1] O. Abidi, M. Heyouni, and K. Jbilou, On some properties of the extended block and global arnoldi methods with applications to model reduction, Numer. Algor., 75 (2014), pp. 285–304.
  * [2] A. Antoulas, Approximation of large-scale dynamical systems, SIAM Philadelphia, (2005).
  * [3] A. Antoulas, C. Beattie, and S. Gugercine, Interpolatory methods for model reduction, Comput. Science and Engineering 21. SIAM Philadelphia, (2020).
  * [4] P. Baranyi, Y. Yam, and P. Várlaki, Tensor product model transformation in polytopic model-based control, (1st ed.). CRC Press, (2014).
  * [5] H. Barkouki, A. H. Bentbib, and K. Jbilou, An extended nonsymmetric block lanczos method for model reduction in large scale dynamical systems, Calcolo, 55(1) (2018), pp. 1–23.
  * [6] A. Barraud, A numerical algorithm to solve atxa - x = q, IEEE Transac. Auto. Contr., 22 (1977), pp. 883–885.
  * [7] P. Benner, M. Ohlberger, A. Cohen, and K. Willcox, Model reduction and approximation, SIAM Philadelphia, (2017).
  * [8] R. Bouyouli, K. Jbilou, R. Sadaka, and H. Sadok, Convergence properties of some block krylov subspace methods for multiple linear systems, J. Comput. Appl. Math., 196 (2006), pp. 498–511.
  * [9] M. J. Brazell, N. LI, C. Navasca, and C. Tamon, Solving multilinear systems via tensor inversion, SIAM J. Matrix Ana. Appl., 34 (2013), pp. 542–570.
  * [10] A. Bunse-Gerstner, D. Kubalińska, G. Vossen, and D. Wilczek, h2-norm optimal model reduction for large scale discrete dynamical mimo systems, J. of Comput. Appl. Math.,, 233(5) (2010), pp. 1202–1216.
  * [11] C. Chen, A. Surana, A. Bloch, and I. Rajapakse, Multilinear time invariant system theory, In 2019 Proceedings of the Conference on Control and its Applications, SIAM, (2019), pp. 118–125.
  * [12] , Data-driven model reduction for multilinear control systems via tensor trains, submitted, https://arxiv.org/pdf/1912.03569.pdf, (2020).
  * [13] , Multilinear time invariant system theory, SIAM J. Contr. Optim.,, 59(1) (2021), pp. 749–776.
  * [14] M. Frangos and I. M. Jaimoukha, Adaptive rational krylov algorithms for model reduction, In Proceedings of the European Control Conference, (2007), pp. 4179–4186.
  * [15] G. H. Golub and C. F. V. Loan, Matrix computations, The Johns Hopkins University Press, (1996).
  * [16] S. Gugercin and A. C. Antoulas, A survey of model reduction by balanced truncation and some new results, Internat. J. Control, 77(8) (2003), pp. 748–766.
  * [17] M. E. Guide, A. E. Ichi, K. Jbilou, and F. Beik, Tensor krylov subspace methods via the einstein product with applications to image and video processing, App. Numer. Math., 181 (2022), pp. 347–363.
  * [18] M. E. Guide, A. E. Ichi, K. Jbilou, and R. Sadaka, On tensor gmres and golub-kahan methods via the t-product for color image processing, The Electronic J. of Linea. Alg., 37 (2021), pp. 524–543.
  * [19] M. A. Hamadi, K. Jbilou, and A. Ratnani, Model reduction method in large scale dynamical systems using an extended-rational block arnoldi method, J. Appl. Math. Comput., 68(1) (2022), pp. 271–293.
  * [20] M. Heyouni, Extended arnoldi methods for large low-rank sylvester matrix equations, Appl. Numer. Math., 60(11) (2010), pp. 1171–1182.
  * [21] M. Heyouni and K. Jbilou, An extended block arnoldi method for large matrix riccati equations, Elect. Trans. Numer. Anal., 33 (2009), pp. 53–62.
  * [22] I.Abdaoui, L. Elbouyahyaoui, and M.Heyouni, An alternative extended block arnoldi method for solving low-rank sylvester equations, Comput. Math. with Appl., 78(8) (2019), pp. 2817–2830.
  * [23] K. Jbilou, A. Messaoudi, and H. Sadok, Global fom and gmres algorithms for matrix equations, Appl. Numer. Math., 31 (1999), pp. 49–63.
  * [24] J. N. Juang and R. S. Pappa, An eigensystem realization algorithm for modal parameter identification and model reduction, J. of Fluid Mechanics, 8 (1985), pp. 620–627.
  * [25] L. Knizhnerman, D. Druskin, and M. Zaslavsky, On optimal convergence rate of the rational krylov subspace reduction for electromagnetic problems in unbounded domains, SIAM J. Numer. Anal., 47(2) (2009), pp. 953–971.
  * [26] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM review, 51 (2009), pp. 455–500.
  * [27] K. Kruppa, Comparison of tensor decomposition methods for simulation of multilinear time-invariant systems with the mti toolbox, IFACPapersOnLine, 50 (2017), pp. 5610–5615.
  * [28] M. Liang and B. Zheng, Further results on moore–penrose inverses of tensors with application to tensor nearness problems, Computers and Math. with App.,, 77(5) (2019), pp. 0898–1221.
  * [29] H. Liu, L. Daniel, and N. Wong, Model reduction and simulation of nonlinear circuits via tensor decomposition., IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst., 34(7) (2015), pp. 1059–1069.
  * [30] V. Mehrmann and T. Stykel, Balanced truncation model reduction for large-scale systems in descriptor form, Lecture Notes in Computational Science and Engineering, 45 (2005), pp. 83–115.
  * [31] B. C. Moore, Principal component analysis in linear systems: controllability, observability and model reduction, IEEE Trans. Auto. Cont., AC-26 (1981), pp. 17–32.
  * [32] M. Nip, J. P. Hespanha, and M. Khammash, Direct numerical solution of algebraic lyapunov equations for large-scale systems using quantized, in 2013 52nd IEEE Conference on Decision and Control (CDC), (2013), pp. 1950–1957.
  * [33] S. Ragnarsson and C. V. Loan, Block tensor unfoldings, SIAM J. Matrix Ana. Appl., 33 (2012), pp. 149–169.
  * [34] M. Rogers, L. Li, and S. J. Russell, Multilinear dynamical systems for tensor time series, in Advances in Neural Information Processing Systems 26 http://papers.nips.cc/paper/5117-multilinear-dynamical-systems-for-tensor-time-series.pdf, 34 (2013), pp. 2634–2642.
  * [35] C. W. Rowley, Model reduction for fluids, using balanced proper orthogonal decomposition, International J. of Bifurcation and Chaos,, 15(3) (2005), pp. 997–1013.
  * [36] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, J. of Fluid Mechanics, 656 (2010), pp. 5–28.
  * [37] V. Simoncini, A new iterative method for solving large-scale lyapunov matrix equations, SIAM J. Sci. Comput.,, 29(3) (2007), pp. 1268–1288.
  * [38] V. Simoncini, D. B. Szyld, and M. Marlliny, On two numerical methods for the solution of large-scale algebraic riccati equations, IMA Journal of Numerical Analysis, 34 (2014), pp. 904–920.
  * [39] A. Surana, G. Patterson, and I. Rajapakse, Dynamic tensor time series modeling and analysis, in 2016 IEEE 55th Conference on Decision and Control (CDC), (2016), pp. 1637–1642.
  * [40] W. Yuchao and W. Yimin, Generalized eigenvalue for even order tensors via einstein product and its applications in multilinear control systems, Comput. Appl. Math.,, 41(8) (2022).