Multilinear Control Systems Theory

Can Chen; Amit Surana; Anthony Bloch; Indika Rajapakse

arXiv:1905.08783·math.OC·March 2, 2021

Multilinear Control Systems Theory

Can Chen, Amit Surana, Anthony Bloch, Indika Rajapakse

PDF

TL;DR

This paper develops a system theoretic framework for multilinear time-invariant systems with tensor states, inputs, and outputs, extending classical LTI concepts using tensor algebra and decompositions for analysis and reduction.

Contribution

It introduces a novel tensor-based approach to MLTI systems, extending stability, reachability, and observability notions, and proposes tensor decomposition methods for model reduction.

Findings

01

Tensor unfolding provides an isomorphism to the general linear group.

02

Tensor decompositions enable efficient system reduction.

03

Numerical examples demonstrate the effectiveness of the proposed methods.

Abstract

In this paper, we provide a system theoretic treatment of a new class of multilinear time-invariant (MLTI) systems in which the states, inputs and outputs are tensors, and the system evolution is governed by multilinear operators. The MLTI system representation is based on the Einstein product and even-order paired tensors. There is a particular tensor unfolding which gives rise to an isomorphism from this tensor space to the general linear group, i.e. the group of invertible matrices. By leveraging this unfolding operation, one can extend classical linear time-invariant (LTI) system notions including stability, reachability and observability to MLTI systems. While the unfolding based formulation is a powerful theoretical construct, the computational advantages of MLTI systems can only be fully realized while working with the tensor form, where hidden patterns/structures can be…

Figures1

Click any figure to enlarge with its caption.

Tables3

Table 1. Table 1: Kronecker rank/TT-ranks approximations of the MLTI system. We omit the first and last trivial TT-ranks in the generalized TTDs of A , B and C .

	Reduced Ranks	# Parameters	$\frac{{‖ G_{full} - G_{red} ‖}_{\infty}}{{‖ G_{full} ‖}_{\infty}}$
Full System	-	783	-
Generalized CPD	$\begin{matrix} 49, 2, 2 \\ 20, 2, 2 \\ 10, 2, 2 \end{matrix}$	$\begin{matrix} 1359 \\ 576 \\ 306 \end{matrix}$	$\begin{matrix} 1.58 \times 10^{- 10} \\ 0.0223 \\ 0.0888 \end{matrix}$
Generalized TTD	$\begin{matrix} {7, 8}, {1, 2}, {2, 2} \\ {7, 6}, {1, 2}, {2, 2} \\ {7, 5}, {1, 2}, {2, 2} \end{matrix}$	$\begin{matrix} 678 \\ 534 \\ 462 \end{matrix}$	$\begin{matrix} 4.39 \times 10^{- 15} \\ 0.0099 \\ 0.4911 \end{matrix}$

Table 2. Table 2: Memory consumption comparison between the generalized TTD and balanced truncation based methods. We reported the TT-ranks of A , B and C (ignoring the first and last trivial TT-ranks) and the number of singular values retained in the Hankel matrix during the balanced truncation.

	Ranks	# Parameters	$\frac{{‖ G_{full} - G_{red} ‖}_{\infty}}{{‖ G_{full} ‖}_{\infty}}$
Full System	-	139968	-
Generalized TTD	{6,6}, {6,6}, {6,6}	5184	$3.98 \times 10^{- 15}$
Balanced Truncation	$\begin{matrix} 200 \\ 100 \\ 40 \end{matrix}$	$\begin{matrix} 120000 \\ 30000 \\ 4800 \end{matrix}$	$\begin{matrix} 0.0169 \\ 0.1001 \\ 0.2360 \end{matrix}$

Table 3. Table 3: Run time comparison between the TTD and SVD based methods in finding the largest singular value of φ ( A ) 𝜑 A \varphi(\textsf{A}) . For the TTD based method, we reported computational time includes conversion from the generalized TTD of A to the TTD of A ~ ~ A \tilde{\textsf{A}} and left- and right-orthonormalization.

$𝒏$	TTD(s)	SVD(s)	$𝝈_{\max}$	Relative error	Stability
6	0.0399	$6.8551 \times 10^{- 4}$	0.8082	$1.3738 \times 10^{- 16}$	asy. stable
8	0.0491	0.0439	0.9626	$4.1523 \times 10^{- 15}$	asy. stable
10	0.0591	0.4979	0.8645	$3.8527 \times 10^{- 15}$	asy. stable
12	0.0909	30.7663	0.8485	$5.7573 \times 10^{- 15}$	asy. stable
14	0.2623	2115.1	0.9984	$1.3566 \times 10^{- 14}$	asy. stable

Equations88

(X \circ Y)_{j_{1} j_{2} \dots j_{N} i_{1} i_{2} \dots i_{M}} = X_{j_{1} j_{2} \dots j_{N}} Y_{i_{1} i_{2} \dots i_{M}} .

(X \circ Y)_{j_{1} j_{2} \dots j_{N} i_{1} i_{2} \dots i_{M}} = X_{j_{1} j_{2} \dots j_{N}} Y_{i_{1} i_{2} \dots i_{M}} .

⟨ X, Y ⟩ = j_{1} = 1 \sum J_{1} \dots j_{N} = 1 \sum J_{N} X_{j_{1} j_{2} \dots j_{N}} Y_{j_{1} j_{2} \dots j_{N}},

⟨ X, Y ⟩ = j_{1} = 1 \sum J_{1} \dots j_{N} = 1 \sum J_{N} X_{j_{1} j_{2} \dots j_{N}} Y_{j_{1} j_{2} \dots j_{N}},

X \times_{1} A_{1} \times_{2} A_{2} \times_{3} \dots \times_{N} A_{N} = X \times {A_{1}, A_{2}, \dots, A_{N}} \in R^{I_{1} \times I_{2} \times \dots \times I_{N}} .

X \times_{1} A_{1} \times_{2} A_{2} \times_{3} \dots \times_{N} A_{N} = X \times {A_{1}, A_{2}, \dots, A_{N}} \in R^{I_{1} \times I_{2} \times \dots \times I_{N}} .

i v ec (j, J) = j_{1} + k = 2 \sum N (j_{k} - 1) l = 1 \prod k - 1 J_{l} .

i v ec (j, J) = j_{1} + k = 2 \sum N (j_{k} - 1) l = 1 \prod k - 1 J_{l} .

(X_{(rc)})_{pq} = X_{p_{1} p_{2} \dots p_{z} q_{1} q_{2} \dots q_{N - z}}^{S},

(X_{(rc)})_{pq} = X_{p_{1} p_{2} \dots p_{z} q_{1} q_{2} \dots q_{N - z}}^{S},

X_{j_{S (1)} j_{S (2)} \dots j_{S (N)}}^{S} = X_{j_{1} j_{2} \dots j_{N}} .

X_{j_{S (1)} j_{S (2)} \dots j_{S (N)}}^{S} = X_{j_{1} j_{2} \dots j_{N}} .

(A * X)_{j_{1} j_{2} \dots j_{N}} = i_{1} = 1 \sum I_{1} \dots i_{N} = 1 \sum I_{N} A_{j_{1} i_{1} \dots j_{N} i_{N}} X_{i_{1} i_{2} \dots i_{N}} .

(A * X)_{j_{1} j_{2} \dots j_{N}} = i_{1} = 1 \sum I_{1} \dots i_{N} = 1 \sum I_{N} A_{j_{1} i_{1} \dots j_{N} i_{N}} X_{i_{1} i_{2} \dots i_{N}} .

(A * B)_{j_{1} i_{1} \dots j_{N} i_{N}} = k_{1} = 1 \sum K_{1} \dots k_{N} = 1 \sum K_{N} A_{j_{1} k_{1} \dots j_{N} k_{N}} B_{k_{1} i_{1} \dots k_{N} i_{N}} .

(A * B)_{j_{1} i_{1} \dots j_{N} i_{N}} = k_{1} = 1 \sum K_{1} \dots k_{N} = 1 \sum K_{N} A_{j_{1} k_{1} \dots j_{N} k_{N}} B_{k_{1} i_{1} \dots k_{N} i_{N}} .

A_{j_{1} i_{1} \dots j_{N} i_{N}} φ A_{i v ec (j, J) i v ec (i, I)},

A_{j_{1} i_{1} \dots j_{N} i_{N}} φ A_{i v ec (j, J) i v ec (i, I)},

{X_{t + 1} = X_{t} \times {A_{1}, \dots, A_{N}} + U_{t} \times {B_{1}, \dots, B_{N}} Y_{t} = X_{t} \times {C_{1}, \dots, C_{N}},

{X_{t + 1} = X_{t} \times {A_{1}, \dots, A_{N}} + U_{t} \times {B_{1}, \dots, B_{N}} Y_{t} = X_{t} \times {C_{1}, \dots, C_{N}},

{X_{t + 1} = A * X_{t} + B * U_{t} Y_{t} = C * X_{t},

{X_{t + 1} = A * X_{t} + B * U_{t} Y_{t} = C * X_{t},

G (z) = C * (z I - A)^{- 1} * B,

G (z) = C * (z I - A)^{- 1} * B,

X_{k} = A^{k} * X_{0} + j = 0 \sum k - 1 A^{k - j - 1} * B * U_{j} .

X_{k} = A^{k} * X_{0} + j = 0 \sum k - 1 A^{k - j - 1} * B * U_{j} .

(A B_{n})_{j_{1} l_{1} \dots j_{N} l_{N}} = ⎩ ⎨ ⎧ A_{j_{1} l_{1} \dots j_{N} l_{N}}, B_{j_{1} l_{1} \dots j_{N} l_{N}}, j_{k} = 1, \dots, J_{k}, l_{k} = 1, \dots, I_{k} \forall k j_{k} = 1, \dots, J_{k} \forall k, l_{k} = 1, \dots, I_{k} for k \neq = n and l_{k} = I_{k} + 1, \dots, 2 I_{k} for k = n .

(A B_{n})_{j_{1} l_{1} \dots j_{N} l_{N}} = ⎩ ⎨ ⎧ A_{j_{1} l_{1} \dots j_{N} l_{N}}, B_{j_{1} l_{1} \dots j_{N} l_{N}}, j_{k} = 1, \dots, J_{k}, l_{k} = 1, \dots, I_{k} \forall k j_{k} = 1, \dots, J_{k} \forall k, l_{k} = 1, \dots, I_{k} for k \neq = n and l_{k} = I_{k} + 1, \dots, 2 I_{k} for k = n .

A B_{n} \in R^{J_{1} \times I_{1} \times \dots \times 2 J_{n} \times I_{n} \times \dots \times J_{N} \times I_{N}}

A B_{n} \in R^{J_{1} \times I_{1} \times \dots \times 2 J_{n} \times I_{n} \times \dots \times J_{N} \times I_{N}}

X = S \times_{1} U_{1} \times_{2} \dots \times_{N} U_{N},

X = S \times_{1} U_{1} \times_{2} \dots \times_{N} U_{N},

X = r = 1 \sum R λ_{r} a_{r}^{(1)} \circ a_{r}^{(2)} \circ \dots \circ a_{r}^{(N)},

X = r = 1 \sum R λ_{r} a_{r}^{(1)} \circ a_{r}^{(2)} \circ \dots \circ a_{r}^{(N)},

X = r_{0} = 1 \sum R_{0} \dots r_{N} = 1 \sum R_{N} X_{r_{0} : r_{1}}^{(1)} \circ X_{r_{1} : r_{2}}^{(2)} \circ \dots \circ X_{r_{N - 1} : r_{N}}^{(N)},

X = r_{0} = 1 \sum R_{0} \dots r_{N} = 1 \sum R_{N} X_{r_{0} : r_{1}}^{(1)} \circ X_{r_{1} : r_{2}}^{(2)} \circ \dots \circ X_{r_{N - 1} : r_{N}}^{(N)},

R_{n} = rank (reshape (X, i = 1 \prod n J_{i}, i = n + 1 \prod N J_{i}))

R_{n} = rank (reshape (X, i = 1 \prod n J_{i}, i = n + 1 \prod N J_{i}))

n = 1 : 2 \sum 2 N k_{A^{(n)}} \geq R + N - 1, n = 2 : 2 \sum 2 N k_{A^{(n)}} \geq R + N - 1

n = 1 : 2 \sum 2 N k_{A^{(n)}} \geq R + N - 1, n = 2 : 2 \sum 2 N k_{A^{(n)}} \geq R + N - 1

W_{r} (t_{0}, t_{1}) = t = t_{0} \sum t_{1} - 1 A^{t_{1} - t - 1} * B * B^{⊤} * (A^{⊤})^{t_{1} - t - 1},

W_{r} (t_{0}, t_{1}) = t = t_{0} \sum t_{1} - 1 A^{t_{1} - t - 1} * B * B^{⊤} * (A^{⊤})^{t_{1} - t - 1},

W_{r} - A * W_{r} * A^{⊤} = B * B^{⊤} .

W_{r} - A * W_{r} * A^{⊤} = B * B^{⊤} .

R = B A * B \dots A^{Π_{J} - 1} * B

R = B A * B \dots A^{Π_{J} - 1} * B

W_{o} (t_{0}, t_{1}) = t = t_{0} \sum t_{1} - 1 (A^{⊤})^{t - t_{0}} * C^{⊤} * C * A^{t - t_{0}},

W_{o} (t_{0}, t_{1}) = t = t_{0} \sum t_{1} - 1 (A^{⊤})^{t - t_{0}} * C^{⊤} * C * A^{t - t_{0}},

A^{⊤} * W_{o} * A - W_{o} = - C^{⊤} * C .

A^{⊤} * W_{o} * A - W_{o} = - C^{⊤} * C .

O = C C * A \dots C * A^{Π_{J} - 1}^{⊤}

O = C C * A \dots C * A^{Π_{J} - 1}^{⊤}

A = r = 1 \sum R A_{r ::}^{(1)} \circ A_{r ::}^{(2)} \circ \dots \circ A_{r ::}^{(N)},

A = r = 1 \sum R A_{r ::}^{(1)} \circ A_{r ::}^{(2)} \circ \dots \circ A_{r ::}^{(N)},

A = r_{0} = 1 \sum R_{0} \dots r_{N} = 1 \sum R_{N} A_{r_{0} :: r_{1}}^{(1)} \circ A_{r_{1} :: r_{2}}^{(2)} \circ \dots \circ A_{r_{N - 1} :: r_{N}}^{(N)},

A = r_{0} = 1 \sum R_{0} \dots r_{N} = 1 \sum R_{N} A_{r_{0} :: r_{1}}^{(1)} \circ A_{r_{1} :: r_{2}}^{(2)} \circ \dots \circ A_{r_{N - 1} :: r_{N}}^{(N)},

A * B = t = 1 \sum T E_{t ::}^{(1)} \circ E_{t ::}^{(2)} \circ \dots \circ E_{t ::}^{(N)},

A * B = t = 1 \sum T E_{t ::}^{(1)} \circ E_{t ::}^{(2)} \circ \dots \circ E_{t ::}^{(N)},

⎩ ⎨ ⎧ X_{t + 1} = r = 1 \sum R_{1} X_{t} \times {A_{r}^{(1)}, \dots, A_{r}^{(N)}} + r = 1 \sum R_{2} U_{t} \times {B_{r}^{(1)}, \dots, B_{r}^{(N)}} Y_{t} = r = 1 \sum R_{3} X_{t} \times {C_{r}^{(1)}, \dots, C_{r}^{(N)}},

⎩ ⎨ ⎧ X_{t + 1} = r = 1 \sum R_{1} X_{t} \times {A_{r}^{(1)}, \dots, A_{r}^{(N)}} + r = 1 \sum R_{2} U_{t} \times {B_{r}^{(1)}, \dots, B_{r}^{(N)}} Y_{t} = r = 1 \sum R_{3} X_{t} \times {C_{r}^{(1)}, \dots, C_{r}^{(N)}},

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

\headers

Multilinear Control Systems TheoryCan Chen, Amit Surana, Anthony Bloch, and Indika Rajapakse

Multilinear Control Systems Theory††thanks: Submitted to the

editors May 17, 2019. This is an extended version of a conference paper presented in SIAM CT19. \fundingThis work is supported in part under AFOSR Award No: FA9550-18-1-0028, NSF grant DMS 1613819, Smale Institute, and the Guaranteeing AI Robustness Against Deception (GARD) program from DARPA (I2O).

Can Chen Department of Mathematics and Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 (). [email protected]

Amit Surana Raytheon Technologies Research Center, East Hartford, CT 06108 (). [email protected]

Anthony Bloch Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 (). [email protected]

Indika Rajapakse Department of Computational Medicine & Bioinformatics, Medical School and Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 (). [email protected]

Abstract

In this paper, we provide a system theoretic treatment of a new class of multilinear time-invariant (MLTI) systems in which the states, inputs and outputs are tensors, and the system evolution is governed by multilinear operators. The MLTI system representation is based on the Einstein product and even-order paired tensors. There is a particular tensor unfolding which gives rise to an isomorphism from this tensor space to the general linear group, i.e. the group of invertible matrices. By leveraging this unfolding operation, one can extend classical linear time-invariant (LTI) system notions including stability, reachability and observability to MLTI systems. While the unfolding based formulation is a powerful theoretical construct, the computational advantages of MLTI systems can only be fully realized while working with the tensor form, where hidden patterns/structures can be exploited for efficient representations and computations. Along these lines, we establish new results which enable one to express tensor unfolding based stability, reachability and observability criteria in terms of more standard notions of tensor ranks/decompositions. In addition, we develop a generalized CANDECOMP/PARAFAC decomposition and tensor train decomposition based model reduction framework, which can significantly reduce the number of MLTI system parameters. We demonstrate our framework with numerical examples.

keywords:

multilinear time-invariant systems, stability, reachability, observability, model reduction, tensor unfolding, tensor ranks/decompositions, block tensors

{AMS}

15A69, 93B05, 93B07, 93B40, 93C05, 93D20

1 Introduction

Controlling high-dimensional systems remains an extremely challenging task as many control strategies do not scale well with the dimension of the systems. Of particular interest in this paper are complex biological and engineering systems in which structure, function and dynamics are highly coupled. Such interactions can be naturally and compactly captured by tensors. Tensors are multidimensional arrays generalized from vectors and matrices, and have wide applications in many domains such as social networks, biology, cognitive science, applied mechanics, scientific computation and signal processing [7, 11, 16, 21, 24]. For example, the organization of the interphase nucleus in the human genome reflects a dynamical interaction between 3D genome structure, function and its relationship to phenotype, a concept known as the 4D Nucleome (4DN) [7]. 4DN research requires a comprehensive view of genome-wide structure, gene expression, the proteome and phenotype, which fits naturally with a tensorial representation [41, 51]. In order to apply the standard system and control framework in applications such as these, tensors need to be vectorized, leading to an extremely high-dimensional system representation in which the number of states/parameters scale exponentially with the number of dimensions of the tensors involved [51]. Moreover, with the vectorization of tensors, hidden patterns/structures, e.g. redundancy/correlations, can get lost, and thus one cannot exploit such inherent structures for efficient representations and computations.

In order to take advantage of tensor algebraic computations, recently a new class of multilinear time-invariant (MLTI) system has been introduced [42, 51], in which the states and outputs are preserved as tensors. The system evolution is generated by the action of multilinear operators which are formed using Tucker products of matrices. By using tensor unfolding, an operation that transforms a tensor into a matrix, Rogers et al. [42] and Surana et al. [51] developed methods for model identification/reduction from tensor time series data. An application of such tensor based representation/identification for skeleton based human behavior recognition from videos demonstrated significant improvements in classification accuracy compared to standard linear time-invariant (LTI) based approaches [11]. However, the MLTI system representation is limited because it assumes the multilinear operators are formed from the Tucker products of matrices (and thus precludes more general tensorial representations) and does not incorporate control inputs.

The role of tensor algebra has also been explored for modeling and simulation of nonlinear dynamics, where the vector field is a multilinear function of states [25]. Tensor decomposition techniques such as CANDECOMP/PARAFAC decomposition (CPD) and tensor train decomposition (TTD) can reduce system size, thus reducing memory usage and enabling efficient computation during simulations. Note that in contrast to the MLTI systems framework of [42, 51], in this application, tensor algebra is applied to the system represented in the conventional vector form. The author in [16] exploits tensor decompositions to compute numerical solutions of master equations associated with Markov processes on extremely large state spaces. The Einstein product and even-order paired tensors, along with TTD, were utilized for developing tensor representations for operators based on nearest-neighbor interactions, construction of pseudoinverses for dimensionality reduction methods and the approximation of transfer operators of dynamical systems.

Similarly, using the Einstein product and even-order paired tensors, Chen et. al. [5] generalized the notion of MLTI systems introduced in [42, 51] and also incorporated control inputs. The Einstein product is a tensor contraction used quite often in tensor calculus and has profound applications in the study of continuum mechanics and the field of relativity theory [14, 27]. Moreover, the space of even-order tensors with the Einstein product has many desirable properties. Brazell et al. [2] discovered that one particular tensor unfolding gives rise to an isomorphism from this tensor space (of even-order tensors equipped with the Einstein product) to the general linear group, i.e. the group of invertible matrices. This isomorphism enables one to define matrix equivalent concepts for tensors including tensor inverse, positive definiteness and eigenvalue decomposition. Using these tensor constructs, Chen et. al. [5] developed tensor algebraic conditions for stability, reachability and observability for generalized input/output MLTI systems. A new notion of block tensors was also introduced which enables one to express these conditions in a compact fashion. Interestingly, these conditions look analogous to the classical conditions for stability, reachability and observability in LTI systems, and reduce to them as a special case.

This paper is an extended version of the introductory paper [5], and in addition to providing various technical details, we also present several new results. The key contributions of this paper are as follows:

In [5], the reachability and observability conditions for MLTI systems were stated in terms of the unfolding rank which requires matricization of the reachability/observability tensors. Here we establish new results relating the unfolding rank to other more standard notions of tensor ranks, including multilinear ranks, CP rank and TT-ranks. Using such relations, we provide criteria for reachability and observability which do not require tensor unfolding, and can be computed using efficient tensor algebraic methods. Similarly, we express MLTI system stability conditions using higher-order singular value decomposition (HOSVD), CPD and TTD. 2. 2.

Using generalized CPD/TTD, we develop a framework for model reduction of MLTI systems. This approach takes advantage of tensor decompositions which otherwise cannot be exploited after unfolding the MLTI systems to obtain a standard LTI form. It also successfully realizes the tensor decomposition based criteria for stability, reachability and observability. Furthermore, we establish new stability results by utilizing the factor matrices from tensor decompositions for this reduced model with lesser computational costs. 3. 3.

We provide computational and memory complexity analysis for the CPD and TTD based methods in comparison to unfolding based matrix methods and demonstrate our framework in four numerical examples.

The paper is organized into nine sections. In Section 2, we review tensor preliminaries including various notions of tensor products, tensor unfolding and properties of even-order paired tensors. Section 3 introduces the MLTI system representation using the Einstein product and even-order paired tensors in detail. In Section 4, we discuss notions of block tensors and tensor ranks/decompositions. We also build new results relating the unfolding rank of a tensor to other more standard notions of tensor ranks. We establish stability, reachability and observability conditions for MLTI systems in Section 5. The application of generalized CPD/TTD for model reduction is discussed in Section 6. Four numerical examples are presented in Section 7. Finally, we summarize different numerical approaches associated with MLTI systems in Section 8 and conclude in Section 9 with future research directions.

2 Tensor preliminaries

We take most of the concepts and notations for tensor algebra from the comprehensive works of Kolda et al. [23, 24] and Ragnarsson et al. [39, 40]. A tensor is a multidimensional array. The order of a tensor is the number of its dimensions, and each dimension is called a mode. An $N$ -th order tensor usually is denoted by $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ . The sets of indexed indices and size of X are denoted by $\textbf{j}=\{j_{1},j_{2},\dots,j_{N}\}$ and $\mathcal{J}=\{J_{1},J_{2},\dots,J_{N}\}$ , respectively. Let $\Pi_{\mathcal{J}}$ represent the product of all elements in $\mathcal{J}$ , i.e. $\Pi_{\mathcal{J}}=\prod_{n=1}^{N}J_{n}$ . It is therefore reasonable to consider scalars $x\in\mathbb{R}$ as zero-order tensors, vectors $\textbf{v}\in\mathbb{R}^{J}$ as first-order tensors, and matrices $\textbf{A}\in\mathbb{R}^{J\times I}$ as second-order tensors.

2.1 Tensor products

By extending the notion of vector outer product, the outer product of two tensors $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ and $\textsf{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{M}}$ is defined as

[TABLE]

In contrast, the inner product of two tensors $\textsf{X},\textsf{Y}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ is defined as

[TABLE]

leading to the tensor Frobenius norm $\|\textsf{X}\|^{2}=\langle\textsf{X},\textsf{X}\rangle$ . We say two tensors X and Y are orthogonal if the inner product $\langle\textsf{X},\textsf{Y}\rangle=0$ . The matrix tensor multiplication $\textsf{X}\times_{n}\textbf{A}$ along mode $n$ for a matrix $\textbf{A}\in\mathbb{R}^{I\times J_{n}}$ is defined by $(\textsf{X}\times_{n}\textbf{A})_{j_{1}j_{2}\dots j_{n-1}ij_{n+1}\dots j_{N}}=\sum_{j_{n}=1}^{J_{n}}\textsf{X}_{j_{1}j_{2}\dots j_{n}\dots j_{N}}\textbf{A}_{ij_{n}}.$ This product can be generalized to what is known as the Tucker product, for $\textbf{A}_{n}\in\mathbb{R}^{I_{n}\times J_{n}}$ ,

[TABLE]

2.2 Tensor unfolding

Tensor unfolding is considered as a critical operation in tensor computations [23, 24, 39]. In order to unfold a tensor $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ into a vector or a matrix, we use an index mapping function $ivec(\cdot,\mathcal{J}):\mathbb{Z}^{+}\times\mathbb{Z}^{+}\times\stackrel{{\scriptstyle\scriptscriptstyle N}}{{\cdots}}\times\mathbb{Z}^{+}\rightarrow\mathbb{Z}^{+}$ as defined by Ragnarsson et al. [39, 40], which is given by

[TABLE]

The index mapping function $ivec$ returns the index for tensor vectorization, i.e. $\textbf{x}\in\mathbb{R}^{\Pi_{\mathcal{J}}}$ is the vectorization of X such that $\textbf{x}_{ivec(\textbf{j},\mathcal{J})}=\textsf{X}_{j_{1}j_{2}\dots j_{N}}$ . If $N=2$ , $ivec$ will stack all the columns of X.

For tensor matricization, let $z$ be an integer such that $1\leq z<N$ , and $\mathbb{S}$ be a permutation of $\{1,2,\dots,N\}$ . If $\textbf{r}=\{\mathbb{S}(1),\mathbb{S}(2),\dots,\mathbb{S}(z)\}$ and $\textbf{c}=\{\mathbb{S}(z+1),\mathbb{S}(z+2),\dots,\mathbb{S}(N)\}$ with $\mathcal{P}=\{J_{\mathbb{S}(1)},J_{\mathbb{S}(2)},\dots,J_{\mathbb{S}(z)}\}$ and $\mathcal{Q}=\{J_{\mathbb{S}(z+1)},J_{\mathbb{S}(z+2)},\dots,J_{\mathbb{S}(N)}\}$ , respectively, the r****c-unfolding matrix of X, denoted by $\textbf{X}_{(\textbf{rc})}\in\mathbb{R}^{\Pi_{\mathcal{P}}\times\Pi_{\mathcal{Q}}}$ , is given by

[TABLE]

where, $p=ivec(\textbf{p},\mathcal{P})$ , $q=ivec(\textbf{q},\mathcal{Q})$ , and $\textsf{X}^{\mathbb{S}}$ is the $\mathbb{S}$ -transpose of X defined as

[TABLE]

When $z=1$ and $\mathbb{S}=\footnotesize{\begin{pmatrix}1&2&\dots&n&n+1&\dots&N\\ n&1&\dots&n-1&n+1&\dots&N\end{pmatrix}}$ , the tensor unfolding is called the $n$ -mode matricization, denoted by $\textbf{X}_{(n)}$ .

2.3 Even-order paired tensors

Here we discuss the notion of even-order paired tensors and the Einstein product which will play an important role in developing the MLTI systems theory.

Definition 2.1.

*Even-order paired tensors are $2N$ -th order tensors with elements specified using a pairwise index notation, i.e. $\textsf{A}_{j_{1}i_{1}\dots j_{N}i_{N}}$ for $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ . *

Definition 2.2.

Given an even-order paired tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , the Einstein product between A and an $N$ -th order tensor $\textsf{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{N}}$ is the contraction along the second index in each pair from A, i.e.

[TABLE]

We use the pairwise index notation for even-order tensors because it is convenient for defining the unfolding transformation $\varphi$ , see Definition 2.4, and for representing core matrices/tensors in tensor decompositions, see Section 6.1. Note that even-order paired tensors and the Einstein product (3) can be viewed as multidimensional generalizations of matrices and the standard matrix-vector product, respectively [16]. Similar to the standard matrix-matrix product, one can also define a generalized form of the Einstein product between two even-order paired tensors. We will see later that the Einstein product can be efficiently computed using tensor decompositions of even-order paired tensors, see Proposition 6.3.

Definition 2.3.

Given two even-order paired tensors $\textsf{A}\in\mathbb{R}^{J_{1}\times K_{1}\times\dots J_{N}\times K_{N}}$ and $\textsf{B}\in\mathbb{R}^{K_{1}\times I_{1}\times\dots\times K_{N}\times I_{N}}$ , the Einstein product $\textsf{A}*\textsf{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ is defined by

[TABLE]

Brazell et al. [2] investigated properties of even-order tensors under the Einstein product (different from (4)) through construction of an isomorphism to GL( $n,\mathbb{R}$ ), i.e. the set of $n\times n$ real valued invertible matrices. The existence of the isomorphism enables one to generalize several matrix concepts, such as invertibility and eigenvalue decomposition to the tensor case [2, 8, 17, 29, 50]. We can establish an analogous isomorphism for even-order paired tensors by a permutation of indices.

Definition 2.4.

Define the map $\varphi$ : $\mathbb{T}_{J_{1}I_{1}\dots J_{N}I_{N}}(\mathbb{R})\rightarrow\mathbb{M}_{\Pi_{\mathcal{J}}\Pi_{\mathcal{I}}}(\mathbb{R})$ with $\varphi(\textsf{A})=\textbf{A}$ defined component-wise as

[TABLE]

*where, $\mathbb{T}_{J_{1}I_{1}\dots J_{N}I_{N}}(\mathbb{R})$ is the set of all real $J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}$ even-order paired tensors, and $\mathbb{M}_{\Pi_{\mathcal{J}}\Pi_{\mathcal{I}}}(\mathbb{R})$ is set of all real $\Pi_{\mathcal{J}}\times\Pi_{\mathcal{I}}$ matrices. *

The map $\varphi$ can be viewed as a tensor unfolding discussed in (2) with $z=N$ and $\mathbb{S}=\footnotesize{\begin{pmatrix}1&2&\dots&N&N+1&N+2&\dots&2N\\ 1&3&\dots&2N-1&2&4&\dots&2N\end{pmatrix}}$ , so the Frobenius norm is preserved through $\varphi$ , i.e. $\|\textsf{A}\|=\|\varphi(\textsf{A})\|$ . More significantly, $\varphi$ is bijective, and the restriction of $\varphi^{-1}$ on the general linear group produces a group isomorphism.

Corollary 2.5.

*Let $J_{n}=I_{n}$ for all $n$ and $\mathbb{G}_{J_{1}J_{1}\dots J_{N}J_{N}}(\mathbb{R})=\varphi^{-1}(\text{GL}(\Pi_{\mathcal{J}},\mathbb{R}))$ , i.e. $\mathbb{G}_{J_{1}J_{1}\dots J_{N}J_{N}}$ is the space of all even-order paired tensors which maps to the general linear group under $\varphi$ . Then $\mathbb{G}_{J_{1}J_{1}\dots J_{N}J_{N}}(\mathbb{R})$ is a group equipped with the Einstein product (4), and $\varphi$ is a group isomorphism. *

Detailed proofs can be found in [2, 17]. Based on the unfolding property, we can define several tensor notations analogous to matrices as follows:

For an even-order paired tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , $\textsf{T}\in\mathbb{R}^{I_{1}\times J_{1}\times\dots\times I_{N}\times J_{N}}$ is called the U-transpose of A if $\textsf{T}_{i_{1}j_{1}\dots i_{N}j_{N}}=\textsf{A}_{j_{1}i_{1}\dots j_{N}i_{N}}$ , and is denoted by $\textsf{A}^{\top}$ . We refer to an even-order paired tensor that is identical to its U-transpose as weakly symmetric. 2. 2.

For an even-order paired tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , the unfolding rank of A is defined as $\text{rank}_{U}(\textsf{A})=\text{rank}(\varphi(\textsf{A}))$ [29]. 3. 3.

An even-order “square” tensor $\textsf{D}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is called the U-diagonal tensor if all its entries are zeros except for $\textsf{D}_{j_{1}j_{1}\dots j_{N}j_{N}}$ . If all the diagonal entires $\textsf{D}_{j_{1}j_{1}\dots j_{N}j_{N}}=1$ , then D is the U-identity tensor, denoted by I. 4. 4.

For an even-order square tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ , if there exists a tensor $\textsf{B}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ such that $\textsf{A}*\textsf{B}=\textsf{B}*\textsf{A}=\textsf{I}$ , then B is called the U-inverse of A, denoted by $\textsf{A}^{-1}$ . 5. 5.

An even-order square tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is called U-positive definite if $\textsf{X}^{\top}*\textsf{A}*\textsf{X}>0$ for any nonzero tensor $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ . 6. 6.

For an even-order square tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ , the unfolding determinant of A is defined as $\text{det}_{U}(\textsf{A})=\text{det}(\varphi(\textsf{A}))$ [29].

In Section A.1, we show that the notion of U-positive definiteness is a generalization of M-positive definiteness and rank-one positive definiteness proposed in [18, 38] for the even-order elasticity tensors.

3 MLTI system representation

To describe the evolution of tensor time series, the authors in [42, 51] introduced a MLTI system using the Tucker product, which can be generalized by incorporating control inputs as follows:

[TABLE]

where, $\textsf{X}_{t}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ is the latent state space tensor, $\textsf{Y}_{t}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{N}}$ is the output tensor, and $\textsf{U}_{t}\in\mathbb{R}^{K_{1}\times K_{2}\times\dots\times K_{N}}$ is a control tensor. $\textbf{A}_{n}\in\mathbb{R}^{J_{n}\times J_{n}}$ , $\textbf{B}_{n}\in\mathbb{R}^{J_{n}\times K_{n}}$ and $\textbf{C}_{n}\in\mathbb{R}^{I_{n}\times J_{n}}$ are real valued matrices for $n=1,2,\dots,N$ . The Tucker product provides a suitable way to deal with MLTI systems because it allows one to exploit matrix computations. However, we find that (6) can be replaced by a more general representation using the notion of even-order paired tensors and the Einstein product. Moreover, the new representation is more concise and systematic compared to the tensor based linear system proposed in [11].

Definition 3.1.

A more general representation of MLTI system is given by

[TABLE]

*where, $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ , $\textsf{B}\in\mathbb{R}^{J_{1}\times K_{1}\times\dots\times J_{N}\times K_{N}}$ and $\textsf{C}\in\mathbb{R}^{I_{1}\times J_{1}\times\dots\times I_{N}\times J_{N}}$ are even-order paired tensors. *

Lemma 3.2.

*Let $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be an even-order paired tensor. Then the product $\textsf{A}\times\{\textbf{U}_{1},\textbf{V}_{1},\dots,\textbf{U}_{N},\textbf{V}_{N}\}=\textsf{U}*\textsf{A}*\textsf{V}^{\top}\in\mathbb{R}^{K_{1}\times L_{1}\times\dots\times K_{N}\times L_{N}}$ for $\textsf{U}=\textbf{U}_{1}\circ\textbf{U}_{2}\circ\dots\circ\textbf{U}_{N}$ and $\textsf{V}=\textbf{V}_{1}\circ\textbf{V}_{2}\circ\dots\circ\textbf{V}_{N}$ where $\textbf{U}_{n}\in\mathbb{R}^{K_{n}\times J_{n}}$ and $\textbf{V}_{n}\in\mathbb{R}^{L_{n}\times I_{n}}$ . *

Proof 3.3.

*This follows from the definitions of the Tucker and Einstein products. *

Proposition 3.4.

*The governing equations (7) can be obtained from (6) by setting A, B and C to be the outer products of component matrices $\{\textbf{A}_{1},\textbf{A}_{2},\dots,\textbf{A}_{N}\}$ , $\{\textbf{B}_{1},\textbf{B}_{2},\dots,\textbf{B}_{N}\}$ and $\{\textbf{C}_{1},\textbf{C}_{2},\dots,\textbf{C}_{N}\}$ , respectively. *

Proof 3.5.

*The result follows from Lemma 3.2 with $I_{n}=1$ and $\textbf{V}_{n}=1$ for all $n$ . *

The main advantages of the MLTI system (7) are as follows:

The Einstein product representation (7) of MLTI systems is indeed the generalization of (6). While Proposition 3.4 shows that MLTI systems in the form of (6) can always be transformed into the form of (7), the converse is not always true, see (23) for example. It is true only when $R_{1}=R_{2}=R_{3}=1$ . 2. 2.

The MLTI system (7) takes a form similar to the standard LTI system model with matrix product replaced with the Einstein product, so the representation is more natural for developing the MLTI systems theory including notions of stability, reachability and observability. Moreover, the concept of transfer functions which is commonly used in modern control theory can be extended for MLTI systems, see Definition 3.6. 3. 3.

We can exploit tensor decompositions (see Section 4.2) of the even-order paired tensors A, B and C to accelerate computations in MLTI systems theory. In particular, if A, B and C possess low tensor rank structures, we can obtain a low-parameter MLTI representation. In addition, many operations such as the Einstein product and unfolding rank can be achieved efficiently in the tensor decomposition format compared to unfolding based matrix methods, see remarks in Sections 5 and 6. 4. 4.

Traditional LTI model reduction and identification techniques such as balanced truncation and eigensystem realization algorithm can be extended using the form of (7).

Definition 3.6.

The transfer function $\textsf{G}(z)$ of (7) is given by

[TABLE]

*where, $z$ is a complex variable. *

We first investigate the elementary solution to the MLTI system (7), which is crucial in the analysis of stability, reachability and observability.

Proposition 3.7.

*For an unforced MLTI system $\textsf{X}_{t+1}=\textsf{A}*\textsf{X}_{t},$ the solution for X at time $k$ , given initial condition $\textsf{X}_{0}$ , is $\textsf{X}_{k}=\textsf{A}^{k}*\textsf{X}_{0}$ where $\textsf{A}^{k}=\textsf{A}*\textsf{A}*\stackrel{{\scriptstyle k}}{{\cdots}}*\textsf{A}$ . *

The proof is straightforward using the notion of even-order paired tensors and the Einstein product. Applying Proposition 3.7, we can write down the explicit solution of (7) which takes an analogous form to the LTI system

[TABLE]

Lastly, we want to note that one can always transform the MLTI system (7) into a LTI system using $\varphi$ , i.e. $\textbf{x}_{t+1}=\varphi(\textsf{A})\textbf{x}_{t}+\varphi(\textsf{B})\textbf{u}_{t}$ , and determine the stability, reachability and observability using classical matrix techniques.

4 Tensor algebra continued

We next discuss notions of block tensors and tensor decompositions which will form the basis for developing tensor algebra based concepts of stability, reachability and observability of the MLTI system (7).

4.1 Block tensors

Analogously to block matrices, one can define the notion of block tensors. For tensors of the same size, we propose a block tensor construction (first appeared in [5]) which does not introduce any wasteful zeros compared to the block tensors proposed in [50], and thus offers computational advantages.

Definition 4.1.

Let $\textsf{A},\textsf{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be two even-order paired tensors. The $n$ -mode row block tensor $\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{n}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{n}\times 2I_{n}\times\dots\times J_{N}\times I_{N}}$ is defined by

[TABLE]

For example, if $\textsf{A},\textsf{B}\in\mathbb{R}^{2\times 2\times 2\times 2}$ , then the 1-mode row block tensor is given by $\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{1}\in\mathbb{R}^{2\times 4\times 2\times 2}$ such that $(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{1})_{:i_{1}::}=\textsf{A}$ for $i_{1}=1,2$ and $(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{1})_{:i_{1}::}=\textsf{B}$ for $i_{1}=3,4$ . Similarly for $\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{2}\in\mathbb{R}^{2\times 2\times 2\times 4}$ . Detailed explanations of the MATLAB colon operation : can be found in Section D.1. When $N=1$ , it reduces to the row block matrices. The $n$ -mode column block tensor

[TABLE]

can be defined in a similar manner. The $n$ -mode block tensors exhibit many properties analogous to block matrix computations, e.g. the Einstein product can distribute over block tensors, and the blocks of $n$ -mode row block tensors map to contiguous blocks under $\varphi$ up to some permutations [39], see details in Section A.2. Therefore, rank is preserved in the block tensor unfolding, i.e. $\text{rank}_{U}(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{n})=\text{rank}(\begin{bmatrix}\varphi(\textsf{A})&\varphi(\textsf{B})\end{bmatrix}),$ where $[\cdot]$ denotes the block matrix operation.

Given $K$ even-order paired tensors $\textsf{X}_{n}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , one can apply Definition 4.1 successively to create a $J_{1}\times I_{1}\times\dots\times J_{n}\times I_{n}K\times\dots\times J_{N}\times I_{N}$ even-order $n$ -mode row block tensor. However, a more general concatenation approach can be defined for multiple blocks.

Definition 4.2.

Given $K$ even-order paired tensors $\textsf{X}_{n}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , if $K=K_{1}K_{2}\dots K_{N}$ , the $J_{1}\times I_{1}K_{1}\times\dots\times J_{N}\times I_{N}K_{N}$ even-order mode row block tensor Y can be constructed in the following way:

Compute the 1-mode row block tensor concatenation over $\{\textsf{X}_{1},\cdots,\textsf{X}_{K_{1}}\}$ , $\{\textsf{X}_{K_{1}+1},\cdots,\textsf{X}_{2K_{1}}\}$ and so on to obtain $K_{2}K_{3}\dots K_{N}$ block tensors denoted by $\textsf{X}_{1}^{(1)},\textsf{X}_{2}^{(1)},\dots,\textsf{X}_{K_{2}K_{3}\dots K_{N}}^{(1)}$ ; 2. 2.

Compute the 2-mode row block tensors concatenation over $\{\textsf{X}_{1}^{(1)},\cdots,\textsf{X}_{K_{2}}^{(1)}\}$ , $\{\textsf{X}_{K_{2}+1}^{(1)},\cdots,\textsf{X}_{2K_{2}}^{(1)}\}$ and so on to obtain $K_{3}K_{4}\dots K_{N}$ block tensors denoted by $\textsf{X}_{1}^{(2)},\textsf{X}_{2}^{(2)},\dots,\textsf{X}_{K_{3}K_{4}\dots K_{N}}^{(2)}$ ; 3. 3.

Keep repeating the process until the last $N$ -mode row block tensor is obtained.

*We denote the mode row block tensor as $\textsf{Y}=\begin{vmatrix}\textsf{X}_{1}&\textsf{X}_{2}&\dots&\textsf{X}_{K}\end{vmatrix}$ . *

For example, suppose that $\textsf{X}_{n}\in\mathbb{R}^{2\times 2\times 2\times 2\times 2\times 2}$ for $n=1,2,\dots,K$ and $K=8$ . Let $K=K_{1}K_{2}K_{3}$ with $K_{1}=K_{2}=K_{3}=2$ . Given this factorization of $K$ , the mode row block tensor $\textsf{Y}\in\mathbb{R}^{2\times 4\times 2\times 4\times 2\times 4}$ is constructed in the manner shown in Fig. 1, in which $\textsf{X}^{(1)}_{n}\in\mathbb{R}^{2\times 4\times 2\times 2\times 2\times 2}$ and $\textsf{X}^{(2)}_{n}\in\mathbb{R}^{2\times 4\times 2\times 4\times 2\times 2}$ . Another factorization with $K_{1}=2$ , $K_{2}=4$ and $K_{3}=1$ would return $\textsf{Y}\in\mathbb{R}^{2\times 4\times 2\times 8\times 2\times 2}$ . The generalized mode column block tensors with multiple blocks can be constructed in a similar manner. When $I_{n}=1$ for all $n$ , the above generalized mode row block tensor maps exactly to contiguous blocks in its unfolding under $\varphi$ , which could be beneficial in many block tensor applications. Furthermore, the choices of $K_{n}$ may affect the structure of mode block tensors, which can be significant in tensor ranks/decompositions [6].

4.2 Tensor ranks and decompositions

There are several definitions of tensor ranks [10, 23, 24], which are intimately related to different notions of tensor decompositions. The multilinear ranks or the $n$ -ranks of X are the ranks of the $n$ -mode matricizations, denoted by $\text{rank}_{n}(\textsf{X})$ . The multilinear ranks are related to the so-called Higher-Order Singular Value Decomposition (HOSVD), a multilinear generalization of the matrix Singular Value Decomposition (SVD) [1, 9].

Theorem 4.3 (HOSVD).

A tensor $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ can be written as

[TABLE]

where, $\textbf{U}_{n}\in\mathbb{R}^{J_{n}\times J_{n}}$ are orthogonal matrices, and $\textsf{S}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ is a tensor of which the subtensors $\textsf{S}_{j_{n}=\alpha}$ obtained by fixing the $n$ -th index to $\alpha$ , have the properties:

all-orthogonality: two subtensors $\textsf{S}_{j_{n}=\alpha}$ and $\textsf{S}_{j_{n}=\beta}$ are orthogonal for all possible values of $n$ , $\alpha$ and $\beta$ subject to $\alpha\neq\beta$ ; 2. 2.

ordering: $\|\textsf{S}_{j_{n}=1}\|\geq\dots\geq\|\textsf{S}_{j_{n}=J_{n}}\|\geq 0$ for all possible values of $n$ .

*The Frobenius norms $\|\textsf{S}_{j_{n}=j}\|$ , denoted by $\gamma_{j}^{(n)}$ , are the $n$ -mode singular values of X. *

De Lathauwer et al. [9] showed that the number of nonvanishing $n$ -mode singular values from the HOSVD of a tensor is equal to its $n$ -mode multilinear rank. In addition, the error bound of the low mutilinear rank approximation is provided in [9]. Unlike the matrix SVD, the approximation fails to obtain the best rank approximation of X. Nevertheless, it still can provide a “good” estimate with appropriate $n$ -mode singular values truncated [9].

Analogous to rank-one matrices, a tensor X is rank-one if it can be written as the outer product of $N$ vectors, i.e. $\textsf{X}=\textbf{a}^{(1)}\circ\textbf{a}^{(2)}\circ\dots\circ\textbf{a}^{(N)}.$ The CANDECOMP/PARAFAC Decomposition (CPD) decomposes a tensor $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ into a sum of rank-one tensors as form of outer products. It is often useful to normalize all the vectors and have weights $\lambda_{r}>0$ in descending order in front:

[TABLE]

where, $\textbf{a}_{r}^{(n)}\in\mathbb{R}^{J_{n}}$ have unit length, and $R$ is called the CP rank of X if it is the minimum integer that achieves (11). The factor matrices $\textbf{A}^{(n)}\in\mathbb{R}^{J_{n}\times R}$ are the combination of the vectors from the rank-one components for $n=1,2,\dots,N$ , i.e. $\textbf{A}^{(n)}=\begin{bmatrix}\textbf{a}_{1}^{(n)}&\textbf{a}_{2}^{(n)}&\dots&\textbf{a}_{R}^{(n)}\end{bmatrix}$ . The CPD is unique up to scaling and permutation under a weak condition: for $N\geq 2$ and $R\geq 2$ , $\sum_{n=1}^{N}k_{\textbf{A}^{(n)}}\geq 2R+(N-1),$ where $k_{\textbf{A}^{(n)}}$ , called the $k$ -rank of a matrix, is the maximum number of columns of $\textbf{A}^{(n)}$ that are linearly independent with each other [26, 47, 49].

The CP rank of a tensor is always greater than or equal to its multilinear ranks [10]. In fact, it is greater than or equal to any unfolding matrix rank [36] (which can be used in unfolding rank and TT-ranks defined later too). The best CP rank approximation is ill-posed [10], but carefully truncating the CP rank will yield a good estimate of the original tensor. Both CPD and HOSVD are special cases of Tucker Decomposition, which decomposes a tensor into the form of Tucker product (1), i.e. $\textsf{Y}=\textsf{X}\times\{\textbf{A}_{1},\textbf{A}_{2},\dots,\textbf{A}_{N}\}$ [25].

The Tensor Train Decomposition (TTD) of $\textsf{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ is given by

[TABLE]

where, $\{R_{0},R_{1},\dots,R_{N}\}$ is the set of TT-ranks with $R_{0}=R_{N}=1$ , and $\textsf{X}^{(n)}\in\mathbb{R}^{R_{n-1}\times J_{n}\times R_{n}}$ are called the core tensors of the TTD [34]. Here we have used : for brevity of notation, for the full definition see Section D.1. Standard TTD algorithms, such as Algorithm 3 in [22], with zero truncation will return the optimal TT-ranks

[TABLE]

for $n=1,2,\dots,N-1$ . A core tensor $\textsf{X}^{(n)}$ is left-orthonormal if $(\bar{\textbf{X}}^{(n)})^{\top}\bar{\textbf{X}}^{(n)}=\textbf{I}\in\mathbb{R}^{R_{n}\times R_{n}}$ , and is right-orthonormal if $\underline{\textbf{X}}^{(n)}(\underline{\textbf{X}}^{(n)})^{\top}=\textbf{I}\in\mathbb{R}^{R_{n-1}\times R_{n-1}}$ , where $\bar{\textbf{X}}^{(n)}=\texttt{reshape}(\textsf{X}^{(n)},R_{n-1}J_{n},R_{n})$ , and $\underline{\textbf{X}}^{(n)}=\texttt{reshape}(\textsf{X}^{(n)},R_{n-1},J_{n}R_{n})$ , respectively [13, 22]. Here reshape refers to the reshape operation in MATLAB, see details in Section D.2. Detailed algorithms for left- and right-orthonormalization can be found in [22]. TTD is advantageous in that it provides better compression, i.e. truncating the TT-ranks results in a quasi-optimal approximation of X, and is computationally more robust [34].

Eigenvalue problems for tensors were first explored by Qi [37] and Lim [28] independently. Brazell et al. [2] formulated a new tensor eigenvalue problem through the isomorphism $\varphi$ for fourth-order tensors, and Cui et al. [8] extended the tensor eigenvalue problem to even-order tensors.

Definition 4.4.

*Let $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ be an even-order square tensor. If $\textsf{X}\in\mathbb{C}^{J_{1}\times J_{2}\dots\times J_{N}}$ is a nonzero $N$ -th order tensor, $\lambda\in\mathbb{C}$ , and X and $\lambda$ satisfy $\textsf{A}*\textsf{X}=\lambda\textsf{X},$ then we call $\lambda$ and X as the U-eigenvalue and U-eigentensor of A, respectively. *

The algebraic and geometric multiplicities of U-eigenvalues can be defined as for matrices. The generalization of the Caley-Hamilton theorem for the tensor case can be obtained by the isomorphism property, i.e. an even-order square tensor A satisfies its own characteristic polynomial $p(\lambda)=\text{det}_{U}(\lambda\textsf{I}-\textsf{A})$ . Moreover, it can be shown that the notion of U-eigenvalues is a generalization of Z-eigenvalues and M-eigenvalues as proposed in [18, 28, 37]. Detailed proofs are omitted in this paper.

Proposition 4.5.

*The tensor eigenvalue problem in Definition 4.4 can be represented by $\textsf{A}=\textsf{V}*\textsf{D}*\textsf{V}^{-1}$ where $\textsf{D}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is an U-diagonal tensor with U-eigenvalues on its diagonal, and $\textsf{V}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is a mode row block tensor consisting of all the U-eigentensors, i.e. $\textsf{V}=\begin{vmatrix}\textsf{X}_{1}&\textsf{X}_{2}&\dots&\textsf{X}_{\Pi_{\mathcal{J}}}\end{vmatrix}.$ We have chosen $K_{n}=J_{n}$ in applying the mode row block tensor operation which enables to express the TEVD in the form analogous to the matrix case. *

Proof 4.6.

*The proof follows immediately from Proposition A.4. *

4.3 Rank relations

We establish new results relating the unfolding rank of an even-order paired tensor to its multilinear ranks, CP rank and TT-ranks. These relationships are useful for checking multilinear generalizations of reachability and observability rank conditions.

Proposition 4.7.

*Let $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be an even-order paired tensor. If $\text{rank}_{U}(\textsf{A})=\Pi_{\mathcal{J}}$ (or $\text{rank}_{U}(\textsf{A})=\Pi_{\mathcal{I}}$ ), then $\text{rank}_{2n-1}(\textsf{A})=J_{n}$ (or $\text{rank}_{2n}(\textsf{A})=I_{n}$ ) for $n=1,2,\dots,N$ . *

Proposition 4.8.

Let $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be an even-order paired tensor given in the CPD format (11) with CP rank equal to $R$ . If the following conditions

[TABLE]

*are satisfied for every $k_{\textbf{A}^{(n)}}\geq 1$ , then $\text{rank}_{U}(\textsf{A})=R$ . *

The notations $\sum_{n=1:2}^{2N}$ and $\sum_{n=2:2}^{2N}$ represent the sums of all odd indices and all even indices, respectively. The detailed proofs of Propositions 4.7 and 4.8 can be found in Sections B.1 and B.2, respectively.

Proposition 4.9.

*Let $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be an even-order paired tensor. Then $\text{rank}_{U}(\textsf{A})=\tilde{R}_{N}$ where $\tilde{R}_{N}$ is the $N$ -th optimal TT-rank of $\tilde{\textsf{A}}$ , the $\mathbb{S}$ -transpose of A with $\mathbb{S}=\footnotesize{\begin{pmatrix}1&2&\dots&N&N+1&N+2&\dots&2N\\ 1&3&\dots&2N-1&2&4&\dots&2N\end{pmatrix}}$ . *

Proof 4.10.

*The result follows from the definition of optimal TT-ranks. *

Remark. Given the TTD of A, the TTD of $\tilde{\textsf{A}}\in\mathbb{R}^{J_{1}\times\dots\times J_{N}\times I_{1}\times\dots\times I_{N}}$ can be constructed by manipulating the core tensors $\textsf{A}^{(n)}$ without converting back to the full format. Assume that $J_{n}=I_{n}=J$ for all $n$ , and $R$ is the average of the TT-ranks of A. If $R$ remains unchanged or decreases during this conversion, the computational complexity is estimated to be at most $\mathcal{O}(N^{2}J^{3}R^{3})$ 111Big O notation: $f(x)=\mathcal{O}(g(x))$ as $x\rightarrow\infty$ if and only if there exists a positive real number $M$ and a real number $x_{0}$ such that $|f(x)|\leq Mg(x)$ for all $x\geq x_{0}$ .. A detailed algorithm for the TTD based permutation can be found in [6].

5 MLTI systems theory

We now introduce the concepts of stability, reachability and observability for MLTI systems. Note that some preliminary results have appeared in our introductory paper [5].

5.1 Stability

There are many notions of stability for dynamical systems [3, 20, 45]. For LTI systems, it is conventional to investigate so-called internal stability. Generalizing from LTI systems, the equilibrium point $\textsf{X}=\textsf{O}$ (O denotes the zero tensors) of an unforced MLTI system is called stable if $\|\textsf{X}_{t}\|\leq\gamma\|\textsf{X}_{0}\|$ for some $\gamma>0$ , asymptotically stable if $\|\textsf{X}_{t}\|\rightarrow 0$ as $t\rightarrow\infty$ , and unstable if it is not stable.

Proposition 5.1.

Let $\lambda_{j}$ be the U-eigenvalues of A for $j=ivec(\textbf{j},\mathcal{J})$ . For an unforced MLTI system, the equilibrium point $\textsf{X}=\textsf{O}$ is:

stable if and only if $|\lambda_{j}|\leq 1$ for all $j=1,2,\dots,\Pi_{\mathcal{J}}$ ; for those equal to 1, its algebraic and geometry multiplicities must be equal; 2. 2.

asymptotically stable if $|\lambda_{j}|<1$ for all $j=1,2,\dots,\Pi_{\mathcal{J}}$ ; 3. 3.

unstable if $|\lambda_{j}|>1$ for some $j=1,2,\dots,\Pi_{\mathcal{J}}$ .

Proof 5.2.

*We only focus on the case when A has a full set of U-eigentensors. It follows from Propositions 4.5 and 3.7 that $\textsf{A}^{k}=\sum_{j_{1}=1}^{J_{1}}\dots\sum_{j_{N}=1}^{J_{N}}\lambda^{k}_{j}\textsf{W}_{j_{1}j_{1}\dots j_{N}j_{N}}$ for some even-order square tensors $\textsf{W}_{j_{1}j_{1}\dots j_{N}j_{N}}$ . Then the results follow immediately. *

Corollary 5.3.

*Suppose that the HOSVD of A is provided with $n$ -mode singular values. For an unforced MLTI system, the equilibrium point $\textsf{X}=\textsf{O}$ is asymptotically stable if the sum of the $n$ -mode singular values square is less than one for any $n$ . *

Proof 5.4.

*Without loss of generality, suppose that $n=1$ . Based on Property 8 in [9], $\sum_{j=1}^{J_{1}}(\gamma_{j}^{(1)})^{2}=\|\textsf{A}\|^{2}=\|\varphi(\textsf{A})\|^{2}.$ In addition, we know that the magnitude of the maximal eigenvalue of a matrix is less than or equal to its Frobenius norm. Hence, the proof follows immediately from Proposition 5.1. *

Corollary 5.5.

*Suppose that the CPD of A is provided and its factor matrices $\textbf{A}^{(n)}$ and $\textbf{A}^{(m)}$ have all the column vectors orthonormal for at least one odd $n$ and even $m$ . For an unforced MLTI system, the equilibrium point $\textsf{X}=\textsf{O}$ is asymptotically stable if the first weight element $\lambda_{1}<1$ . *

The proof of Corollary 5.5 is presented in Section B.3.

Corollary 5.6.

*Suppose that the TTD of $\tilde{\textsf{A}}\in\mathbb{R}^{J_{1}\times\dots\times J_{N}\times J_{1}\times\dots\times J_{N}}$ , defined in Proposition 4.9, is provided with the first $N-1$ core tensors left-orthonormal and the last $N$ core tensors right-orthonormal. For an unforced MLTI system, the equilibrium point $\textsf{X}=\textsf{O}$ is asymptotically stable if the largest singular value of $\bar{\tilde{\textbf{A}}}^{(N)}$ is less than one, where $\bar{\tilde{\textbf{A}}}^{(N)}=\texttt{reshape}(\tilde{\textsf{A}}^{(N)},R_{N-1}J_{N},R_{N})$ . *

Proof 5.7.

*Based on the results of [22], the singular values of $\bar{\tilde{\textbf{A}}}^{(N)}$ are the singular values of $\varphi(\textsf{A})$ . In addition, we know that the magnitude of the maximal eigenvalue of a matrix is less than or equal to its largest singular value. Hence, the proof follows immediately from Proposition 5.1. *

Remark. Although Proposition 5.1 offers strong stability results for unforced MLTI systems, computing U-eigenvalues usually requires an order of $\mathcal{O}(\Pi_{\mathcal{J}}^{3})$ number of operations through tensor unfolding and matrix eigenvalue decomposition. To the contrary, Corollaries 5.3, 5.5 and 5.6 can be used to determine the stability of MLTI systems much faster. In particular, if the TTD of $\tilde{\textsf{A}}$ is provided, the time complexity of left- and right-orthonormalization is about $\mathcal{O}(NJR^{3})$ assuming that $J_{n}=J$ for all $n$ , and $R$ is the average of the TT-ranks of $\tilde{\textsf{A}}$ [34]. Moreover, truncating the TT-rank $\tilde{R}_{N}$ of $\tilde{\textsf{A}}$ would not alter the largest singular values of $\bar{\tilde{\textbf{A}}}^{(N)}$ . Therefore, setting $\tilde{R}_{N}=1$ and computing the vector 2-norm of $\bar{\tilde{\textbf{A}}}^{(N)}$ will return the largest singular value of $\varphi(\textsf{A})$ .

5.2 Reachability

In this and the following subsections, we introduce the definitions of reachability and observability for MLTI systems which are similar to analogous concepts for the LTI systems [3, 20, 45]. We then establish sufficient and necessary conditions for reachability and observability for MLTI systems.

Definition 5.8.

*The MLTI system (7) is said to be reachable on $[t_{0},t_{1}]$ if, given any initial condition $\textsf{X}_{0}$ and any final state $\textsf{X}_{1}$ , there exists a sequence of inputs $\textsf{U}_{t}$ that steers the state of the system from $\textsf{X}_{t_{0}}=\textsf{X}_{0}$ to $\textsf{X}_{t_{1}}=\textsf{X}_{1}$ . *

Theorem 5.9.

The pair $(\textsf{A},\textsf{B})$ is reachable on $[t_{0},t_{1}]$ if and only if the reachability Gramian

[TABLE]

*which is a weakly symmetric even-order square tensor, is U-positive definite. *

Proof 5.10.

Suppose $\textsf{W}_{r}(t_{0},t_{1})$ is U-positive definite, and let $\textsf{X}_{0}$ be the initial state and $\textsf{X}_{1}$ be the desired final state. Choose $\textsf{U}_{t}=\textsf{B}^{\top}*(\textsf{A}^{\top})^{t_{1}-t-1}*\textsf{W}_{r}^{-1}(t_{0},t_{1})*\textsf{V}$ for some constant tensor V. It follows from the solution of the system (7) that $\textsf{X}_{t_{1}}=\textsf{A}^{t_{1}}*\textsf{X}_{0}+\sum_{j=0}^{t_{1}-1}\textsf{A}^{t_{1}-j-1}*\textsf{B}*\textsf{U}_{t}=\textsf{A}^{t_{1}}*\textsf{X}_{0}+\textsf{W}_{r}(t_{0},t_{1})*\textsf{W}_{r}^{-1}(t_{0},t_{1})*\textsf{V}=\textsf{A}^{t_{1}}*\textsf{X}_{0}+\textsf{V}.$ Taking $\textsf{V}=-\textsf{A}^{t_{1}}*\textsf{X}_{0}+\textsf{X}_{1}$ , we have $\textsf{X}_{t_{1}}=\textsf{X}_{1}$ .

*We show the converse by contradiction. Suppose $\textsf{W}_{r}(t_{0},t_{1})$ is not U-positive definite. Then there exists $\textsf{X}_{a}\neq\textsf{O}$ such that $\textsf{X}_{a}^{\top}*\textsf{A}^{t_{1}-t-1}*\textsf{B}=\textsf{O}$ for any $t$ . Take $\textsf{X}_{1}=\textsf{X}_{a}+\textsf{A}^{t_{1}}*\textsf{X}_{0}$ , and it follows that $\textsf{X}_{a}+\textsf{A}^{t_{1}}*\textsf{X}_{0}=\textsf{A}^{t_{1}}*\textsf{X}_{0}+\sum_{j=t_{0}}^{t_{1}-1}\textsf{A}^{t_{1}-j-1}*\textsf{B}*\textsf{U}_{j}\,.$ Multiplying from the left by $\textsf{X}_{a}^{\top}$ yields $\textsf{X}_{a}^{\top}*\textsf{X}_{a}=\sum_{j=t_{0}}^{t_{1}-1}\textsf{X}_{a}^{\top}*\textsf{A}^{t_{1}-j-1}*\textsf{B}*\textsf{U}_{j}=0,$ which implies that $\textsf{X}_{a}=\textsf{O}$ , a contradiction. *

Corollary 5.11.

*If the reachability Gramian $\textsf{W}_{r}(t_{0},t_{1})$ is not M-positive definite, the pair $(\textsf{A},\textsf{B})$ is not reachable on $[t_{0},t_{1}]$ . *

Proof 5.12.

*The proof follows immediately from Propositions A.2 and 5.9. *

The reachability Gramian assesses to what degree each state is affected by an input [44]. The infinite horizon reachability Gramian can be computed from the tensor Lyapunov equation which is defined by

[TABLE]

By the unfolding property, if the pair $(\textsf{A},\textsf{B})$ is reachable over an infinite horizon and all the U-eigenvalues of A have magnitude less than one, one can show that there exists a unique weakly symmetric U-positive definite solution $\textsf{W}_{r}$ . Solving the infinite horizon reachability Gramian from the tensor Lyapunov equation may be computationally intensive, so a tensor version of the Kalman rank condition is also provided.

Proposition 5.13.

The pair $(\textsf{A},\textsf{B})$ is reachable if and only if the $J_{1}\times J_{1}K_{1}\times\dots\times J_{N}\times J_{N}K_{N}$ even-order reachability tensor

[TABLE]

*spans $\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ . In other words, $\text{rank}_{U}(\mathscr{R})=\Pi_{\mathcal{J}}$ . *

Proof 5.14.

*The proof follows from Proposition A.4 and the generalized Cayley Hamilton theorem discussed in the tensor eigenvalue problem. *

First, any choice of construction for the mode row block tensor works for the reachability tensor. Second, when $N=1$ , Proposition 5.13 simplifies to the famous Kalman rank condition for reachability of LTI systems. The following corollaries involving with HOSVD (multilinear ranks), CPD (CP rank) and TTD (TT-ranks) provide useful necessary or sufficient conditions for reachability of MLTI systems if the reachability tensor $\mathscr{R}$ is given in the HOSVD, CPD or TTD format.

Corollary 5.15.

*Given the reachability tensor $\mathscr{R}$ in (16), if $\text{rank}_{2n-1}(\mathscr{R})\neq J_{n}$ for some $n$ , the pair $(\textsf{A},\textsf{B})$ is not reachable. *

Proof 5.16.

*The proof follows immediately from Propositions 4.7 and 5.13. *

Corollary 5.17.

*Given the reachability tensor $\mathscr{R}$ in (16), if the set of $n$ -mode singular values of $\mathscr{R}$ obtained from the HOSVD contains zero for odd $n$ , the pair $(\textsf{A},\textsf{B})$ is not reachable. *

Proof 5.18.

*We know that the number of nonvanishing $n$ -mode singular values equals to its corresponding $n$ -mode multilinear rank. Hence, the result follows immediately from Propositions 5.13 and 5.15. *

Corollary 5.19.

*Given the reachability tensor $\mathscr{R}$ in (16), if the CPD of $\mathscr{R}$ satisfies (13) with CP rank equal to $\Pi_{\mathcal{J}}$ , the pair $(\textsf{A},\textsf{B})$ is reachable. Conversely, if the pair $(\textsf{A},\textsf{B})$ is reachable, then the CP rank of $\mathscr{R}$ is greater than or equal to $\Pi_{\mathcal{J}}$ . *

Proof 5.20.

*The first part of the proof follows immediately from Propositions 4.8 and 5.13. The second part of the proof follows from the fact that the CP rank of a tensor is greater than or equal to its unfolding rank. *

Corollary 5.21.

*Given the reachability tensor $\mathscr{R}$ in (16), the pair $(\textsf{A},\textsf{B})$ is reachable if and only if the $N$ -th optimal TT-rank of $\tilde{\mathscr{R}}\in\mathbb{R}^{J_{1}\times\dots\times J_{N}\times J_{1}K_{1}\times\dots\times J_{N}K_{N}}$ , defined in Proposition 4.9, is equal to $\Pi_{\mathcal{J}}$ . *

Proof 5.22.

*The proof follows immediately from Propositions 4.9 and 5.13. *

Remark. Finding the unfolding rank of the reachability tensor $\mathscr{R}$ through tensor unfolding and matrix QR decomposition is computationally expensive, and has a $\mathcal{O}(\Pi_{\mathcal{J}}^{3}\Pi_{\mathcal{K}})$ time complexity. However, if the reachability tensor $\mathscr{R}$ is already given in the tensor decomposition format, computing the unfolding rank can be achieved efficiently based on Corollaries 5.17, 5.19 and 5.21. Particularly, if the TTD of $\tilde{\mathscr{R}}$ is provided, we do not need any additional computation to obtain the unfolding rank.

5.3 Observability

The results of observability can be simply obtained by the duality principle, similarly to LTI systems.

Definition 5.23.

*The MLTI system (6) is said to be observable on $[t_{0},t_{1}]$ if any initial state $\textsf{X}_{t_{0}}=\textsf{X}_{0}$ can be uniquely determined by $\textsf{Y}_{t}$ on $[t_{0},t_{1}]$ . *

Theorem 5.24.

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

[TABLE]

*which is a weakly symmetric even-order square tensor, is U-positive definite. *

Proof 5.25.

Suppose that $\textsf{W}_{o}(t_{0},t_{1})$ is U-positive definite and let $\textsf{X}_{0}$ be the initial state such that $\textsf{Y}_{t}=\textsf{C}*\textsf{X}_{t}=\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{0}$ for any $t\in[t_{0},t_{1}]$ . Multiplying from the left by $(\textsf{A}^{\top})^{t-t_{0}}*\textsf{C}^{\top}$ yields $(\textsf{A}^{\top})^{t-t_{0}}*\textsf{C}^{\top}*\textsf{Y}_{t}=(\textsf{A}^{\top})^{t-t_{0}}*\textsf{C}^{\top}*\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{0}$ , which implies that $\sum_{t=t_{0}}^{t_{1}-1}(\textsf{A}^{\top})^{t-t_{0}}*\textsf{C}^{\top}*\textsf{Y}_{t}=\sum_{t=t_{0}}^{t_{1}-1}(\textsf{A}^{\top})^{t-t_{0}}*\textsf{C}^{\top}*\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{0}=\textsf{W}_{o}(t_{0},t_{1})*\textsf{X}_{0}.$ Since $\textsf{W}_{o}(t_{0},t_{1})$ is U-invertible, this equation has a unique solution $\textsf{X}_{0}=\textsf{W}_{o}^{-1}(t_{0},t_{1})\sum_{t=t_{0}}^{t_{1}-1}(\textsf{A}^{\top})^{t-t_{0}}*\textsf{C}^{\top}*\textsf{Y}_{t}.$ Hence, $(\textsf{A},\textsf{C})$ is observable on $[t_{0},t_{1}]$ .

*Again, we show the converse by contradiction. Suppose that $\textsf{W}_{o}(t_{0},t_{1})$ is not U-positive definite. Then there exists $\textsf{X}_{a}\neq\textsf{O}$ such that $\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{a}=\textsf{O}$ for any $t$ . Take $\textsf{X}_{t_{0}}=\textsf{X}_{0}+\textsf{X}_{a}$ for some initial state $\textsf{X}_{0}$ . Then $\textsf{Y}_{t}=\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{0}+\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{a}=\textsf{C}*\textsf{A}^{t-t_{0}}*\textsf{X}_{0}$ for any $t\in[t_{0},t_{1}]$ . The initial states $\textsf{X}_{0}$ and $\textsf{X}_{0}+\textsf{X}_{a}$ produce the same output, which implies that $(\textsf{A},\textsf{C})$ is not observable on $[t_{0},t_{1}]$ , a contradiction. *

Corollary 5.26.

*If the observability Gramian $\textsf{W}_{o}(t_{0},t_{1})$ is not M-positive definite, the pair $(\textsf{A},\textsf{C})$ is not observable on $[t_{0},t_{1}]$ . *

The observability Gramian assesses to what degree each state affects future outputs [44]. The infinite horizon observability Gramian can be computed from the tensor Lyapunov equation defined by

[TABLE]

If the pair $(\textsf{A},\textsf{C})$ is observable and all the U-eigenvalues of A have magnitude less than one, there exists a unique weakly symmetric U-positive definite solution $\textsf{W}_{o}$ .

The following results can be proved similarly to those in Section 5.2.

Proposition 5.27.

The pair $(\textsf{A},\textsf{C})$ is observable if and only if the $I_{1}J_{1}\times J_{1}\times\dots\times I_{N}J_{N}\times J_{N}$ even-order observability tensor

[TABLE]

*spans $\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ . In other words, $\text{rank}_{U}(\mathscr{O})=\Pi_{\mathcal{J}}$ . *

Corollary 5.28.

*Given the observability tensor $\mathscr{O}$ in (19), if $\text{rank}_{2n}(\mathscr{O})\neq J_{n}$ for some $n$ , the pair $(\textsf{A},\textsf{C})$ is not observable. *

Corollary 5.29.

*Given the observability tensor $\mathscr{O}$ in (19), if the set of $n$ -mode singular values of $\mathscr{O}$ obtained from the HOSVD contains zero for even $n$ , the pair $(\textsf{A},\textsf{C})$ is not observable. *

Corollary 5.30.

*Given the observability tensor $\mathscr{O}$ in (19), if the CPD of $\mathscr{O}$ satisfies (13) with CP rank equal to $\Pi_{\mathcal{J}}$ , the pair $(\textsf{A},\textsf{C})$ is observable. Conversely, if the pair $(\textsf{A},\textsf{C})$ is observable, then the CP rank of $\mathscr{O}$ is greater than or equal to $\Pi_{\mathcal{J}}$ . *

Corollary 5.31.

*Given the observability tensor $\mathscr{O}$ in (19), the pair $(\textsf{A},\textsf{C})$ is observable if and only if the $N$ -th optimal TT-rank of $\tilde{\mathscr{O}}\in\mathbb{R}^{I_{1}J_{1}\times\dots\times I_{N}J_{N}\times J_{1}\times\dots\times J_{N}}$ , defined in Proposition 4.9, is equal to $\Pi_{\mathcal{J}}$ . *

6 Model reduction for MLTI systems

Based on the observations in Section 5, it is more natural to manipulate MLTI systems in the tensor decomposition format so that all the computational advantages can be realized. This may also result in a more compressed representation.

6.1 Generalized CPD/TTD

We first introduce the notion of generalized CPD/TTD for even-order paired tensors described in [16], in which the generalized CPD can also be viewed as the extension of the Kronecker rank approximation proposed by Van Loan [52]. Generalized CPD and TTD share a similar format and possess many analogous properties.

Definition 6.1.

Given an even-order paired tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , the generalized CPD of A is given by

[TABLE]

*where, $\textsf{A}^{(n)}\in\mathbb{R}^{R\times J_{n}\times I_{n}}$ . Extending Van Loan’s definition [52], we call the smallest $R$ that achieves (20) the Kronecker rank of A. *

Definition 6.2.

Given an even-order paired tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ , the generalized TTD of A is given by

[TABLE]

*where, $\textsf{A}^{(n)}\in\mathbb{R}^{R_{n-1}\times J_{n}\times I_{n}\times R_{n}}$ , and $\{R_{0},R_{1},\dots,R_{N}\}$ is the set of TT-ranks with $R_{0}=R_{N}=1$ . *

Please refer to Section D.1 for the use of : notation. Given two even-order paired tensors in the generalized CPD/TTD format, the Einstein product (4) between the two can be computed without having to reconstruct the full tensors, i.e. keeping the original format [16]. The following proposition states the case for generalized CPD, which also applies to generalized TTD.

Proposition 6.3.

Given two even-order paired tensors $\textsf{A}\in\mathbb{R}^{J_{1}\times K_{1}\times\dots J_{N}\times K_{N}}$ and $\textsf{B}\in\mathbb{R}^{K_{1}\times I_{1}\times\dots\times K_{N}\times I_{N}}$ in the format of (20) with Kronekcer ranks $R$ and $S$ , respectively, the Einstein product $\textsf{A}*\textsf{B}$ is given by

[TABLE]

*where, $\textsf{E}^{(n)}_{t::}=\textsf{A}^{(n)}_{r::}\textsf{B}^{(n)}_{s::}\in\mathbb{R}^{J_{n}\times I_{n}}$ , and $t=ivec(\{r,s\},\{R,S\})$ with $T=RS$ . *

Remark. The computational complexity of the Einstein product (22) is about $\mathcal{O}(NJ^{3}R^{2})$ assuming that $J_{n}=I_{n}=K_{n}=J$ and $R=S$ , which is much lower than $\mathcal{O}(J^{3N})$ from the Einstein product (4) if $R$ is small.

The generalized CPD can be recovered from the standard CPD, and similarly for generalized TTD (see Algorithm 2). The algorithm below is extended from the results by Van Loan [52] about the Kronecker rank approximation. Thus, one can easily obtain generalized CPD by using any technique for computing the standard CPD including alternating least square (ALS) and modified ALS methods [23, 24].

6.2 MLTI model reduction

The problem of model reduction has been studied heavily in the framework of classical control [12, 15, 33]. Methods including proper orthogonal decomposition (POD), scale-separation and averaging, balanced truncation are applied in many engineering applications when dealing with high-dimensional linear/nonlinear systems [32]. As mentioned in Section 3, using generalized CPD/TTD, we propose a new MLTI representation with fewer parameters. Note that we omit colons in each component tensor in this and the following subsections for simplicity (e.g. $\textsf{A}_{r}^{(n)}=\textsf{A}_{r::}^{(n)}$ ).

Proposition 6.4.

The MLTI system (7) is equivalent to

[TABLE]

*where, $R_{1},R_{2},R_{3}$ are the Kronecker ranks of the system, and $\textsf{A}^{(n)}\in\mathbb{R}^{R_{1}\times J_{n}\times J_{n}}$ , $\textsf{B}^{(n)}\in\mathbb{R}^{R_{2}\times J_{n}\times K_{n}}$ and $\textsf{C}^{(n)}\in\mathbb{R}^{R_{3}\times I_{n}\times J_{n}}$ . *

Proof 6.5.

*The proof follows from Definitions 6.1 and 3.4. *

Remark. The number of parameters of the MLTI system representation (23) is $R_{1}\sum_{n=1}^{N}J_{n}^{2}+R_{2}\sum_{n=1}^{N}J_{n}K_{n}+R_{3}\sum_{n=1}^{N}I_{n}J_{n}$ . If the Kronecker ranks $R_{1},R_{2},R_{3}$ are relatively small, the total number of parameters is much less than that of the MLTI system model (7) which is given by $\prod_{n=1}^{N}J_{n}^{2}+\prod_{n=1}^{N}J_{n}K_{n}+\prod_{n=1}^{N}I_{n}J_{n}$ .

The MLTI system representation (23) is attractive for systems captured by sparse tensors or tensors with low Kronecker ranks where the two advantages, model reduction and computational efficiency, can be exploited. In particular, if A, B and C are fourth-order paired tensors, the generalized CPDs are reduced to matrix SVD problems, see section 9.2 in [52]. However, there are two major drawbacks. First, for $N>2$ , there is no exact method to compute the Kronecker rank of a tensor [23], and truncating the rank does not ensure a good estimate. Second, current CPD algorithms are not numerically stable, which could result in ill-conditioning during the tensor decomposition and low rank approximation. One way to fix these issues is to replace generalized CPD by generalized TTD in (23), which takes a similar form. The algorithms for computing generalized TTD are numerically stable with unique optimal TT-ranks [34]. Most importantly, the TTD based results obtained in Section 5 can be realized in the form of (23). For example, we can determine the stability of MLTI systems from the TTD of $\tilde{\textsf{A}}$ defined in Proposition 4.9, which can be obtained from the generalized TTD of A efficiently (similar to the remark in Section 4.3).

Recall from Section 3 that one can always convert the MLTI system (7) to an equivalent LTI form and then apply traditional model reduction approaches, e.g. balanced truncation. However, after converting to a matrix form, the low tensor rank structure exploited in the form of (23) may not be preserved, and thus low memory requirements cannot be achieved, see Section 7.3. Furthermore, as shown in [6], the MLTI system (23) can be used to further develop a higher-order balanced truncation framework directly in the TTD format, which can provide additional computation and memory benefits over unfolding based model reduction methods.

6.3 Explicit solution and stability

In addition to using tensor decompositions, we can exploit matrix calculations of the factor matrices $\textsf{A}^{(n)}_{r}$ to develop notions including explicit solution and stability for the MLTI system (23) which also have lower computational costs compared to unfolding based methods.

Proposition 6.6 (Solution).

For an unforced MLTI system $\textsf{X}_{t+1}=\sum_{r=1}^{R_{1}}\textsf{X}_{t}\times\{\textsf{A}_{r}^{(1)},\textsf{A}_{r}^{(2)},\dots,\textsf{A}_{r}^{(N)}\},$ the solution for X at time $k$ , given initial condition $\textsf{X}_{0}$ , is

[TABLE]

*where, $\bar{\textsf{A}}_{r}^{(n)}=\textsf{A}_{r_{1}}^{(n)}\textsf{A}_{r_{2}}^{(n)}\dots\textsf{A}_{r_{k}}^{(n)}$ for $r=ivec(\{r_{1},r_{2},\dots,r_{k}\},\{R_{1},R_{1},\stackrel{{\scriptstyle k}}{{\cdots}},R_{1}\})$ . *

Proof 6.7.

*The result follows immediately from Propositions 6.3 and 3.4. *

If the Kronecker rank $R_{1}$ is small, computing the explicit solution using (24) can be faster than using the Einstein product (4). Additionally, we can assess the stability of the unforced MLTI system of (23) based upon the Lyapunov approach.

Proposition 6.8 (Stability).

For the unforced MLTI system of (23), the equilibrium point $\textsf{X}=\textsf{O}$ is

stable (i.s.L) if $\sum_{r=1}^{R_{1}}\prod_{n=1}^{N}\alpha_{r}^{(n)}=1$ ; 2. 2.

asymptotically stable (i.s.L) if $\sum_{r=1}^{R_{1}}\prod_{n=1}^{N}\alpha_{r}^{(n)}<1$ ,

*where, $\alpha_{r}^{(n)}$ denote the largest singular values of $\textsf{A}_{r}^{(n)}$ . *

Proof 6.9.

*Let’s consider $V(\textsf{X})=\|\textsf{X}\|$ as the Lyapunov function candidate and let $f(\textsf{X})=\sum_{r=1}^{R_{1}}\textsf{X}\times\{\textsf{A}_{r}^{(1)},\textsf{A}_{r}^{(2)},\dots,\textsf{A}_{r}^{(N)}\}$ . Then it follows that $V(f(\textsf{X}))-V(\textsf{X})=\|\sum_{r=1}^{R_{1}}\textsf{X}\times\{\textsf{A}_{r}^{(1)},\textsf{A}_{r}^{(2)},\dots,\textsf{A}_{r}^{(N)}\}\|-\|\textsf{X}\|\leq\sum_{r=1}^{R_{1}}\|\textsf{X}\times\{\textsf{A}_{r}^{(1)},\textsf{A}_{r}^{(2)},\dots,\textsf{A}_{r}^{(N)}\}\|-\|\textsf{X}\|\leq(\sum_{r=1}^{R_{1}}\prod_{n=1}^{N}\alpha_{r}^{(n)}-1)\|\textsf{X}\|,$ where the last inequality is based on Theorem 6 in [19]. Then the results follow immediately. *

Remark. The computational complexity of finding the matrix SVDs of the factor matrices can be estimated as $\mathcal{O}(NJ^{3}R_{1})$ assuming that $J_{n}=J$ for all $n$ .

When all the Kronecker ranks of the system $R_{1}=R_{2}=R_{3}=1$ , the MLTI system (23) reduces to the Tucker product representation proposed by Surana et al. [51], which provides a more direct way to see that the Tucker based MLTI model is only a special case of the MLTI system (7). Additionally, we can obtain stronger stability conditions for the unforced MLTI system in this case.

Proposition 6.10 (Stability).

Suppose that $R_{1}=1$ in (23), and $\rho^{(n)}$ are the spectral radii of $\textsf{A}_{1}^{(n)}$ . Then the unforced MLTI system of (23) is

stable if and only if $\prod_{n=1}^{N}\rho^{(n)}\leq 1$ , and when $\prod_{n=1}^{N}\rho^{(n)}=1$ , their corresponding eigenvalues must have equal algebraic and geometric multiplicity; 2. 2.

asymptotically stable if $\prod_{n=1}^{N}\rho^{(n)}<1$ ; 3. 3.

unstable if $\prod_{n=1}^{N}\rho^{(n)}>1$ .

Proof 6.11.

*Based on Equation (2.25) in [39], $\varphi(\textsf{A})=\textsf{A}_{1}^{(N)}\otimes\textsf{A}_{1}^{(N-1)}\otimes\dots\otimes\textsf{A}_{1}^{(1)}$ where the operation $\otimes$ denotes the Kronecker product. Moreover, the U-eigenvalues of A are equal to the products of eigenvalues of these component matrices $\textsf{A}_{1}^{(n)}$ , and the U-eigenvalues have equal algebraic and geometric multiplicities if and only if the factor eigenvalues have equal multiplicities [4]. Then the results follow immediately from Proposition 5.1. *

The above results including Propositions 6.4, 6.6 and 6.8 can be reformulated by replacing the Kronecker rank summation by a series of TT-ranks summations if A, B and C are given in the generalized TTD format. Finally, the Kronecker product can be used to unfold the MLTI system (23) into a LTI system, i.e.

[TABLE]

and similarly for $\varphi(\textsf{B})$ and $\varphi(\textsf{C})$ . Hence, one can apply traditional control theory techniques to determine the MLTI system properties.

7 Numerical examples

We provide four examples to illustrate the MLTI systems theory and model reduction using the techniques developed above. All the numerical examples presented were performed on a Linux machine with 8 GB RAM and a 2.4 GHz Intel Core i5 processor and were conducted in MATLAB R2018a with the Tensor Toolbox 2.6 [46] and the TT toolbox [35].

7.1 Reachability and observability tensors

In this example, we consider a simple single-input and single-output (SISO) system that is given by (6) with $\textbf{A}_{1}=\begin{bmatrix}0&1&0\\ 0&0&1\\ 0.2&0.5&0.8\end{bmatrix},\text{ }\textbf{A}_{2}=\begin{bmatrix}0&1\\ 0.5&0\end{bmatrix},\textbf{B}_{1}=\begin{bmatrix}0\\ 0\\ 1\end{bmatrix},\text{ }\textbf{B}_{2}=\begin{bmatrix}0\\ 1\end{bmatrix},\textbf{C}_{1}=\begin{bmatrix}1&0&0\end{bmatrix},\text{ }\textbf{C}_{2}=\begin{bmatrix}1&0\end{bmatrix},$ and the states $\textsf{X}_{t}\in\mathbb{R}^{3\times 2}$ are second-order tensors, i.e. matrices. The product of the two spectral radii of $\textbf{A}_{1}$ and $\textbf{A}_{2}$ is 0.9207, which implies that the system is asymptotically stable. In addition, the reachability and observability tensors based on (16) and (19) are given by

[TABLE]

and

[TABLE]

respectively. We compute the TTDs of the permuted tensors $\tilde{\mathscr{R}}$ and $\tilde{\mathscr{O}}$ , respectively and observe that $\text{rank}_{U}(\mathscr{R})=6$ and $\text{rank}_{U}(\mathscr{O})=6$ . The system therefore is both reachable and observable.

7.2 Kronecker rank/TT-ranks approximation

In this example, we consider a SISO MLTI system (7) with random sparse tensors $\textsf{A}\in\mathbb{R}^{3\times 3\times 3\times 3\times 3\times 3}$ , $\textsf{B}\in\mathbb{R}^{3\times 3\times 3}$ and $\textsf{C}\in\mathbb{R}^{3\times 3\times 3}$ . According to Algorithm 1, we compute the generalized CPDs of A, B and C using the tensor toolbox function cp_als with estimated Kronecker ranks $R_{1}=49$ , $R_{2}=2$ and $R_{3}=2$ , respectively, see Generalized CPD in Table 1. Note that the number of parameters in the system with full Kronekcer ranks could be greater than that for the original system. We then fix $R_{2}$ and $R_{3}$ and gradually truncate $R_{1}$ , since $R_{1}$ is most critical in determining the number of parameters of the reduced system. As we can see in the table, the number of parameters decreases dramatically as $R_{1}$ decreases. In order to assess the approximation error resulting from this truncation, we compute the relative error using the $\mathcal{H}$ -infinity norm $\|\cdot\|_{\infty}$ between the full system and reduced system transfer functions based on (8). In particular, we find that when $R_{1}=10$ , the reduced MLTI system is still close to the original system with $\mathcal{H}$ -infinity norm relative error of 0.0888.

We repeat a similar process for TT-ranks approximation through generalized TTD (see Algorithm 2). The results are shown in the same table. We find that both generalized CPD and TTD can achieve efficient model reduction while keeping the approximation errors low. Generalized TTD in particular achieves better accuracy for a similar number of reduced parameters as compared to generalized CPD, but the latter can maintain a resonable approximation error with an even lower number of parameters. The Bode diagrams for the reduced MLTI systems are shown in Fig. 2. Note that in this example, we manually selected the truncation to study the tradeoff between number of parameters in the reduced system and the approximation error.

7.3 Memory consumption comparison

In this example, we consider a multiple-input and multiple-output (MIMO) MLTI system (7) with random even-order paired tensors $\textsf{A},\textsf{B},\textsf{C}\in\mathbb{R}^{6\times 6\times 6\times 6\times 6\times 6}$ that possess low TT-ranks. We compare the memory consumptions of the generalized TTD based representation (23) with the reduced models obtained from the unfolding based balanced truncation. The results are shown in Table 2. One can clearly see that if the MLTI systems possess low TT-ranks structure, the generalized TTD based approach achieves much better accuracy for a similar number of parameters as compared to balanced truncation.

7.4 Computational time comparison

In this example, we consider unforced MLTI systems (23) with random sparse even-order paired tensors $\textsf{A}\in\mathbb{R}^{2\times 2\times\stackrel{{\scriptstyle\scriptscriptstyle n}}{{\cdots}}\times 2\times 2}$ in the generalized TTD format such that $\varphi(\textsf{A})\in\mathbb{R}^{2^{n}\times 2^{n}}$ . We compare the run time of Corollary 5.6 with the matrix SVD of $\varphi(\textsf{A})$ for determining the stability of the systems. The results are shown in Table 3. When $n\geq 10$ , the TTD based method for finding the largest singular value of $\varphi(\textsf{A})$ exhibits a signficant time advantage compared to the matrix SVD based method for which the time increases exponentially.

8 Discussion

While tensor unfolding to a matrix form provides the advantage of leveraging highly optimized matrix algebra libraries, in doing so however one may not be able to exploit the higher-order hidden patterns/structures, e.g. redundancy/correlations, present in the tensor. For instance, in the context of solving PDEs, Brazell et al. [2] found that higher-order tensor representations preserve low bandwidth, thereby keeping the computational cost and memory requirement low. As shown in Sections 7.4 and 7.3, TTD based methods are more efficient in terms of computational speed and memory requirements compared to unfolding based methods when the MLTI systems have low TT-ranks structure. Although CPD typically offers better compression than TTD, the computation of CP rank is NP-hard, and the lower rank approximations can be ill-posed [10]. TTD is more suitable for numerical computations with well developed TT-algebra [34]. Basic tensor operations such as addition, the Einstein product, Frobenius norm, block tensor, solution to multilinear equations and tensor pseudoinverse, can be computed and maintained in the TTD format, without requiring full tensor representation. This can provide significant computational advantages in finding the reachability/observability tensors and associated unfolding ranks according to Corollaries 5.21 and 5.31, and in obtaining solution of the tensor Lyapunov equations. For details, we refer the reader to [6] and the references therein.

Another line of approach is to exploit the isomorphism property to build algorithms directly in the full tensor format from existing methods. For example, Brazell et al. [2] proposed Higher-Order Biconjugate Gradient (HOBG) method for solving multilinear systems which can be used for computing U-inverses and MLTI system transfer functions. Analogously, one can generalize the matrix based Rayleigh Quotient Iteration method for computing U-eigenvalues (which can be used for determining MLTI system stability) directly in the tensor form, see Algorithm 3. However, the computational efficiency of this type of method remains to be investigated. Finally, one can combine tensor algebra based and matrix based methods to provide the advantages of both approaches as hybrid methods, see some examples in [6] in the context of MLTI model reduction. In future, it would be worthwhile to systematically explore which of the above mentioned approaches or combination thereof is best given the problem structure.

9 Conclusion

In this paper, we provided a comprehensive treatment of a newly introduced MLTI system representation using even-order paired tensors and the Einstein product. We established new results which enable one to express tensor unfolding based stability, reachability and observability criteria in terms of more standard notions of tensor ranks/decompositions. We introduced a generalized CPD/TTD based model reduction framework which can significantly reduce the number of MLTI system parameters and realize the tensor decomposition based methods. We also presented computational complexity analysis of our proposed framework, and illustrated the benefits through numerical examples. In particular, TTD offers several computational advantages over CPD and HOSVD, and provides a good representational choice for facilitating numerical computations associated with MLTI systems.

As mentioned in Section 8, more work is required to fully realize the potential of tensor algebra based computations for MLTI systems. It will also be worthwhile to develop theoretical and computational framework for observer and feedback control design for MLTI systems, and apply these techniques in real world complex systems. One particular application we plan to investigate is that of cellular reprogramming which involves introducing transcription factors as a control mechanism to transform one cell type to another. These systems naturally have matrix or tensor state spaces describing their genome-wide structure and gene expression [30, 43]. Such applications would also need to account for nonlinearity and stochasticity in tensor based dynamical system representation and analysis framework, and is an important direction for future research.

Appendix A Additional tensor algebra

A.1 M-positive definiteness/rank-one positive definiteness

Definition A.1.

An even-order square tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is called M-positive definite if the multilinear functional

[TABLE]

*for any nonzero vector $\textbf{x}_{n}$ . If all $\textbf{x}_{n}$ are equal, A is called rank-one positive definite. *

Proposition A.2.

*If an even-order square tensor $\textsf{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is U-positive definite, it is M-positive definite. Moreover, if $J_{1}=J_{2}=\dots=J_{N}$ , U-positive definiteness also implies rank-one positive definiteness. *

Proof A.3.

*By Lemma 3.2, it follows that $\textsf{A}\times\{\textbf{x}_{1}^{\top},\textbf{x}_{1}^{\top},\dots,\textbf{x}_{N}^{\top},\textbf{x}_{N}^{\top}\}=\textsf{X}^{\top}*\textsf{A}*\textsf{X}$ for $\textsf{X}=\textbf{x}_{1}\circ\textbf{x}_{2}\circ\dots\circ\textbf{x}_{N}$ , i.e. X is a rank-one tensor. Moreover, if $J_{1}=J_{2}=\dots=J_{N}$ , M-positive definiteness implies rank-one positive definiteness [38]. Therefore, the results follow immediately. *

A.2 Block tensor properties

Proposition A.4.

Let $\textsf{A},\textsf{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ and $\textsf{C},\textsf{D}\in\mathbb{R}^{I_{1}\times K_{1}\times\dots\times I_{N}\times K_{N}}$ . Then the following properties hold:

$\mathcal{P}*\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{n}=\begin{vmatrix}\textsf{P}*\textsf{A}&\textsf{P}*\textsf{B}\end{vmatrix}_{n}$ * for any $\textsf{P}\in\mathbb{R}^{L_{1}\times J_{1}\times\dots\times L_{N}\times J_{N}}$ ;* 2. 2.

$\begin{vmatrix}\textsf{C}\\ \textsf{D}\end{vmatrix}_{n}*\textsf{Q}=\begin{vmatrix}\textsf{C}*\textsf{Q}\\ \textsf{D}*\textsf{Q}\end{vmatrix}_{n}$ * for any $\textsf{Q}\in\mathbb{R}^{K_{1}\times R_{1}\times\dots\times K_{N}\times R_{N}}$ ;* 3. 3.

$\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{n}*\begin{vmatrix}\textsf{C}\\ \textsf{D}\end{vmatrix}_{n}=\textsf{A}*\textsf{C}+\textsf{B}*\textsf{D}$ .

Proof A.5.

*The proof follows immediately from the definition of $n$ -mode row/column block tensors and the Einstein product. *

Proposition A.6.

*Let $\textsf{A},\textsf{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be two even-order paired tensors. Then $\varphi(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{n})=\begin{bmatrix}\varphi(\textsf{A})&\varphi(\textsf{B})\end{bmatrix}\textbf{P},$ where P is a column permutation matrix. In particular, when $I_{n}=1$ for all $n$ or $n=N$ , P is the identity matrix. *

Proof A.7.

We consider the case for $N=2$ . Since the size of the odd modes of the block tensor remains the same, we only need to consider the even modes’ unfolding transformation. When $n=1$ , the index mapping function for the even modes is

[TABLE]

for $i_{1}=1,2,\dots,2I_{1}$ . Based on the definition of $n$ -mode row block tensors, the first $I_{1}$ columns of $\varphi(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{1})$ are the vectorizations of $\textsf{A}_{:i_{1}:i_{2}}$ for $i_{1}=1,2,\dots,I_{1}$ and $i_{2}=1$ , and the second $I_{1}$ columns are the vectorizations of $\textsf{B}_{:i_{1}:i_{2}}$ for $i_{1}=I_{1}+1,I_{1}+2,\dots,2I_{1}$ and $i_{2}=1$ . The alternating pattern continues for all $I_{2}$ pairs of $I_{1}$ columns. Hence, $\varphi(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{1})=\begin{bmatrix}\varphi(\textsf{A})&\varphi(\textsf{B})\end{bmatrix}\textbf{P}$ for some column permutation matrix P. When $n=2$ , the index mapping function for the even modes is given by

[TABLE]

*for $i_{2}=1,2,\dots,2I_{2}$ . Similarly, the first $I_{1}I_{2}$ columns of $\varphi(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{2})$ are the vectorizations of $\textsf{A}_{:i_{1}:i_{2}}$ for $i_{1}=1,2,\dots,I_{1}$ and $i_{2}=1,2,\dots,I_{2}$ , and the second $I_{1}I_{2}$ columns are the vectorizations of $\textsf{B}_{:i_{1}:i_{2}}$ for $i_{1}=1,2,\dots,I_{1}$ and $i_{2}=I_{2}+1,I_{2}+2,\dots,2I_{2}$ . Hence, $\varphi(\begin{vmatrix}\textsf{A}&\textsf{B}\end{vmatrix}_{2})=\begin{bmatrix}\varphi(\textsf{A})&\varphi(\textsf{B})\end{bmatrix}$ . A similar analysis can be used to prove the case for $N>2$ . Moreover, when $I_{n}=1$ for all $n$ , $\varphi(\textsf{A})$ and $\varphi(\textsf{B})$ are vectors, so no permutation needs to be considered. The proposition can be considered as a special case of Theorem 3.3 in [39]. *

Appendix B Tensor ranks/decompositions proofs

B.1 Proof of Proposition 4.7

Without loss of generality, assume that $\Pi_{\mathcal{I}}\leq\Pi_{\mathcal{J}}$ and $\text{rank}_{U}(\textsf{A})=\Pi_{\mathcal{I}}$ . Then $\varphi(\textsf{A})$ has $\Pi_{\mathcal{I}}$ linearly independent columns. The goal here is to construct a transformation from $\varphi(\textsf{A})$ to $\textbf{A}_{(2n)}^{\top}$ , which can be easily visualized through the representation $(z,\mathbb{S})$ defined in (2). Let

[TABLE]

Clearly, $\varphi(\textsf{A})$ and $\textbf{A}_{(2n)}^{\top}$ can be represented by $(N,\mathbb{S}_{1})$ and $(2N-1,\mathbb{S}_{2})$ , respectively. According to the definition of the index mapping function $ivec(\textbf{i},\mathcal{I})$ , we first require a column permutation matrix P such that $\varphi(\textsf{A})\textbf{P}$ is represented by $(N,\mathbb{S}_{3})$ . Every $I_{n}$ columns of $\varphi(\textsf{A})\textbf{P}$ correspond to the columns of $\textbf{A}_{(2n)}^{\top}$ . Collect each set of $I_{n}$ columns of $\varphi(\textsf{A})\textbf{P}$ and stack them vertically to form a tall matrix $\tilde{\textbf{A}}$ with the representation $(2N-1,\mathbb{S}_{4})$ . Since the columns of $\varphi(\textsf{A})\textbf{P}$ are linearly independent, $\text{rank}(\tilde{\textbf{A}})=I_{n}$ . Finally, according to the definition of the index mapping function $ivec(\textbf{j},\mathcal{J})$ , we require a row permutation matrix Q such that $\textbf{Q}\tilde{\textbf{A}}=\textbf{A}_{(2n)}^{\top}$ . Hence, $\text{rank}_{2n}(\textsf{A})=\text{rank}(\textbf{A}_{(2n)}^{\top})=I_{n}$ . Note that the converse of the statement is incorrect.

B.2 Proof of Proposition 4.8

In order to prove Proposition 4.8, we need to introduce the concept of Khatri-Rao product.

Definition B.1.

Given two matrices $\textbf{A}\in\mathbb{R}^{J\times I}$ and $\textbf{B}\in\mathbb{R}^{K\times I}$ , the Khatri-Rao product, denoted by $\textbf{A}\odot\textbf{B}$ , results in a $JK\times I$ matrix:

[TABLE]

*where, $\otimes$ denotes the Kronecker product, and $\textbf{a}_{n}$ and $\textbf{b}_{n}$ are the column vectors of A and B, respectively. *

The following lemma provided by Sidiropoulos et al. [48, 49] gives some properties of rank and $k$ -rank of the Khatri-Rao product $\textbf{A}\odot\textbf{B}$ .

Lemma B.2.

*Given two matrices $\textbf{A}\in\mathbb{R}^{J\times R},\textbf{B}\in\mathbb{R}^{I\times R}$ , the Khatri-Rao product $\textbf{A}\odot\textbf{B}$ has column rank $R$ if $k_{\textbf{A}}+k_{\textbf{B}}\geq R+1$ for $k_{\textbf{A}},k_{\textbf{B}}\geq 1$ . Moreover, $k_{\textbf{A}\odot\textbf{B}}\geq\min{\{k_{\textbf{A}}+k_{\textbf{B}}-1,R\}}$ . *

Proposition B.3.

*Given matrices $\textbf{A}^{(n)}\in\mathbb{R}^{J_{n}\times R}$ , the Khatri-Rao product $\textbf{A}^{(1)}\odot\textbf{A}^{(2)}\odot\dots\odot\textbf{A}^{(N)}$ has column rank $R$ if $\sum_{n=1}^{N}k_{\textbf{A}^{(n)}}\geq R+N-1$ for $k_{\textbf{A}^{(n)}}\geq 1$ . *

Proof B.4.

Suppose that $N=3$ . By Lemma B.2, the Khatri-Rao product $\textbf{A}^{(1)}\odot\textbf{A}^{(2)}\odot\textbf{A}^{(3)}$ has full column rank $R$ if $k_{\textbf{A}^{(1)}\odot\textbf{A}^{(2)}}+k_{\textbf{A}^{(3)}}\geq R+1.$ Since we know that $k_{\textbf{A}\odot\textbf{B}}\geq\min{\{k_{\textbf{A}}+k_{\textbf{B}}-1,R\}}$ , the above inequality can be satisfied if

[TABLE]

*When $k_{\textbf{A}^{(1)}}+k_{\textbf{A}^{(2)}}>R+1$ , the condition is reduced to $k_{\textbf{A}^{(3)}}\geq 1$ , and when $k_{\textbf{A}^{(1)}}+k_{\textbf{A}^{(2)}}\leq R+1$ , the condition becomes $k_{\textbf{A}^{(1)}}+k_{\textbf{A}^{(2)}}+k_{\textbf{A}^{(3)}}\geq R+2$ . Therefore, the Khatri-Rao product $\textbf{A}^{(1)}\odot\textbf{A}^{(2)}\odot\textbf{A}^{(3)}$ has full column rank $R$ if $k_{\textbf{A}^{(1)}}+k_{\textbf{A}^{(2)}}+k_{\textbf{A}^{(3)}}\geq R+2.$ The result can be easily extended to $n=N$ using the same approach. *

Now, we can prove Proposition 4.8. Suppose that A has the CPD format (11) with CP rank equal to $R$ . Applying the unfolding transformation $\varphi$ yields

[TABLE]

where, $\textbf{S}\in\mathbb{R}^{R\times R}$ is a diagonal matrix containing the weights of the CPD on its diagonal. By Proposition B.3, the two Khatri-Rao products $\textbf{A}^{(2N-1)}\odot\dots\odot\textbf{A}^{(1)}$ and $\textbf{A}^{(2N)}\odot\dots\odot\textbf{A}^{(2)}$ have full column rank $R$ if the two conditions $\sum_{n=1:2}^{2N}k_{\textbf{A}^{(n)}}\geq R+N-1\text{, and }\sum_{n=2:2}^{2N}k_{\textbf{A}^{(n)}}\geq R+N-1$ are satisfied. Hence, $\text{rank}_{U}(\textsf{A})=R$ . Note that we do not require the CPD of A to be unique in the statement.

B.3 Proof of Corollary 5.5

The proof is formulated similarly to the one above. We need to use the properties of Khatri-Rao product.

Lemma B.5.

*Given matrices $\textbf{A}^{(n)}\in\mathbb{R}^{J_{n}\times R}$ , the Khatri-Rao product $\textbf{A}^{(1)}\odot\textbf{A}^{(2)}\odot\dots\odot\textbf{A}^{(N)}$ has all the column vectors orthogonal if at least one of $\textbf{A}^{(n)}$ has all the column vectors orthogonal for $n=1,2,\dots,N$ . *

Proof B.6.

Suppose that $N=2$ . Based on the properties of Kronecker product, for any $1\leq n,m\leq R$ , the inner product between $\textbf{a}_{n}^{(1)}\otimes\textbf{a}_{n}^{(2)}$ and $\textbf{a}_{m}^{(1)}\otimes\textbf{a}_{m}^{(2)}$ is given by

[TABLE]

*Therefore, if $\textbf{A}^{(1)}$ or $\textbf{A}^{(2)}$ has all column vectors orthogonal, then the inner product between $\textbf{a}_{n}^{(1)}\otimes\textbf{a}_{n}^{(2)}$ and $\textbf{a}_{m}^{(1)}\otimes\textbf{a}_{m}^{(2)}$ is zero for any $n,m$ . *

Now we can prove Corollary 5.5. Suppose that A has the CPD format (11). Applying the unfolding transformation $\varphi$ yields

[TABLE]

where, $\textbf{S}\in\mathbb{R}^{R\times R}$ is a diagonal matrix containing the weights of the CPD on its diagonal. By Lemma B.5, the two Khatri-Rao products $\textbf{A}^{(2N-1)}\odot\dots\odot\textbf{A}^{(1)}$ and $\textbf{A}^{(2N)}\odot\dots\odot\textbf{A}^{(2)}$ have all the column vectors orthonormal if $\textbf{A}^{(n)}$ and $\textbf{A}^{(m)}$ have all the column vectors orthonormal for at least one odd $n$ and even $m$ . Thus, $\lambda_{1}$ will be the largest singular value of $\varphi(\textsf{A})$ . In addition, we know that the magnitude of the maximal eigenvalue of a matrix is less than or equal to its largest singular value. Hence, the proof follows immediately from Proposition 5.1. Note that there is one special case when the CPD uniqueness condition fails, i.e. $\sum_{n=1}^{2N}k_{\textbf{A}^{(n)}}=2R+2N-2$ . However, different CPDs, satisfying the orthonormal condition, correspond to the same matrix SVD under $\varphi$ up to some orthogonal transformations.

Appendix C Numerical algorithms

Appendix D MATLAB functions

D.1 The colon operator

The colon : is one of the most useful operators in MATLAB, which can create vectors, subscript arrays and specify for iterations. For our purpose, it acts as shorthand to include all subscripts in a particular array dimension [31]. For example, $\textbf{A}_{:i}$ is equivalent to $\textbf{A}_{ji}$ for all $j$ . In the following, we represent TTD, generalized CPD and TTD in the component-wise form.

(12) $\Leftrightarrow\textsf{X}_{j_{1}j_{2}\dots j_{N}}=\sum_{r_{0}=1}^{R_{0}}\dots\sum_{r_{N}=1}^{R_{N}}\textsf{X}^{(1)}_{r_{0}j_{1}r_{1}}\textsf{X}^{(2)}_{r_{1}j_{2}r_{2}}\dots\textsf{X}^{(N)}_{r_{N-1}j_{N}r_{N}}$ . 2. 2.

(20) $\Leftrightarrow\textsf{A}_{j_{1}i_{1}\dots j_{N}i_{N}}=\sum_{r=1}^{R}\textsf{A}_{rj_{1}i_{1}}^{(1)}\textsf{A}_{rj_{2}i_{2}}^{(2)}\dots\textsf{A}_{rj_{N}i_{N}}^{(N)}$ . 3. 3.

(21) $\Leftrightarrow\textsf{A}_{j_{1}i_{1}\dots j_{N}i_{N}}=\sum_{r_{0}=1}^{R_{0}}\dots\sum_{r_{N}=1}^{R_{N}}\textsf{A}^{(1)}_{r_{0}j_{1}i_{1}r_{1}}\textsf{A}^{(2)}_{r_{1}j_{2}i_{2}r_{2}}\dots\textsf{A}^{(N)}_{r_{N-1}j_{N}i_{N}r_{N}}$ .

D.2 The reshape operator

The command $\textsf{B}=\texttt{reshape}(\textsf{A},J_{1},J_{2},\dots,J_{N})$ reshapes a tensor A into a $J_{1}\times J_{2}\times\dots\times J_{N}$ order tensor such that the number of elements in B matches the number of elements in A [31].

Acknowledgments

We would like to thank Dr. Frederick Leve at the Air Force Office of Scientific Research (AFOSR) for support and encouragement. We would also like to thank the two referees for their constructive comments, which lead to a significant improvement of the paper.

Bibliography52

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] G. Bergqvist and E. G. Larsson , The higher-order singular value decomposition: Theory and an application [lecture notes] , IEEE Signal Processing Magazine, 27 (2010), pp. 151–154, https://doi.org/10.1109/MSP.2010.936030 . · doi ↗
2[2] M. J. Brazell, N. Li, C. Navasca, and C. Tamon , Solving multilinear systems via tensor inversion , SIAM Journal on Matrix Analysis Applications, 34 (2013), pp. 542–570, https://doi.org/10.1137/100804577 . · doi ↗
3[3] R. Brockett , Finite Dimensional Linear Systems , Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015, https://doi.org/10.1137/1.9781611973884 . · doi ↗
4[4] B. J. Broxson , The kronecker product , UNF Theses and Dissertations, (2006), http://digitalcommons.unf.edu/etd/25 .
5[5] C. Chen, A. Surana, A. Bloch, and I. Rajapakse , Multilinear time invariant system theory , in 2019 Proceedings of SIAM Conference on Control and its Applications, pp. 118–125, https://doi.org/10.1137/1.9781611975758.18 . · doi ↗
6[6] C. Chen, A. Surana, A. Bloch, and I. Rajapakse , Data-driven model reduction for multilinear control systems via tensor trains , 2019, https://arxiv.org/abs/1912.03569 .
7[7] H. Chen, J. Chen, L. A. Muir, S. Ronquist, W. Meixner, M. Ljungman, T. Ried, S. Smale, and I. Rajapakse , Functional organization of the human 4d nucleome , Proceedings of the National Academy of Sciences, 112 (2015), pp. 8002–8007, https://doi.org/10.1073/pnas.1505822112 . · doi ↗
8[8] L.-B. Cui, C. Chen, W. Li, and M. K. Ng , An eigenvalue problem for even order tensors with its applications , Linear and Multilinear Algebra, 64 (2016), pp. 602–621, https://doi.org/10.1080/03081087.2015.1071311 . · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Multilinear Control Systems Theory††thanks: Submitted to the

Abstract

keywords:

1 Introduction

2 Tensor preliminaries

2.1 Tensor products

2.2 Tensor unfolding

2.3 Even-order paired tensors

Definition 2.1**.**

Definition 2.2**.**

Definition 2.3**.**

Definition 2.4**.**

Corollary 2.5**.**

3 MLTI system representation

Definition 3.1**.**

Lemma 3.2**.**

Proof 3.3**.**

Proposition 3.4**.**

Proof 3.5**.**

Definition 3.6**.**

Proposition 3.7**.**

4 Tensor algebra continued

4.1 Block tensors

Definition 4.1**.**

Definition 4.2**.**

4.2 Tensor ranks and decompositions

Theorem 4.3** (HOSVD).**

Definition 4.4**.**

Proposition 4.5**.**

Proof 4.6**.**

4.3 Rank relations

Proposition 4.7**.**

Proposition 4.8**.**

Proposition 4.9**.**

Proof 4.10**.**

5 MLTI systems theory

5.1 Stability

Proposition 5.1**.**

Proof 5.2**.**

Corollary 5.3**.**

Proof 5.4**.**

Corollary 5.5**.**

Corollary 5.6**.**

Proof 5.7**.**

5.2 Reachability

Definition 5.8**.**

Theorem 5.9**.**

Proof 5.10**.**

Corollary 5.11**.**

Proof 5.12**.**

Proposition 5.13**.**

Proof 5.14**.**

Corollary 5.15**.**

Proof 5.16**.**

Corollary 5.17**.**

Proof 5.18**.**

Corollary 5.19**.**

Proof 5.20**.**

Corollary 5.21**.**

Proof 5.22**.**

5.3 Observability

Definition 5.23**.**

Theorem 5.24**.**

Proof 5.25**.**

Corollary 5.26**.**

Proposition 5.27**.**

Corollary 5.28**.**

Corollary 5.29**.**

Corollary 5.30**.**

Corollary 5.31**.**

6 Model reduction for MLTI systems

6.1 Generalized CPD/TTD

Definition 6.1**.**

Definition 6.2**.**

Definition 2.1.

Definition 2.2.

Definition 2.3.

Definition 2.4.

Corollary 2.5.

Definition 3.1.

Lemma 3.2.

Proof 3.3.

Proposition 3.4.

Proof 3.5.

Definition 3.6.

Proposition 3.7.

Definition 4.1.

Definition 4.2.

Theorem 4.3 (HOSVD).

Definition 4.4.

Proposition 4.5.

Proof 4.6.

Proposition 4.7.

Proposition 4.8.

Proposition 4.9.

Proof 4.10.

Proposition 5.1.

Proof 5.2.

Corollary 5.3.

Proof 5.4.

Corollary 5.5.

Corollary 5.6.

Proof 5.7.

Definition 5.8.

Theorem 5.9.

Proof 5.10.

Corollary 5.11.

Proof 5.12.

Proposition 5.13.

Proof 5.14.

Corollary 5.15.

Proof 5.16.

Corollary 5.17.

Proof 5.18.

Corollary 5.19.

Proof 5.20.

Corollary 5.21.

Proof 5.22.

Definition 5.23.

Theorem 5.24.

Proof 5.25.

Corollary 5.26.

Proposition 5.27.

Corollary 5.28.

Corollary 5.29.

Corollary 5.30.

Corollary 5.31.

Definition 6.1.

Definition 6.2.

Proposition 6.3.

Proposition 6.4.

Proof 6.5.

Proposition 6.6 (Solution).

Proof 6.7.

Proposition 6.8 (Stability).

Proof 6.9.

Proposition 6.10 (Stability).

Proof 6.11.

Definition A.1.

Proposition A.2.

Proof A.3.

Proposition A.4.

Proof A.5.

Proposition A.6.

Proof A.7.

Definition B.1.

Lemma B.2.

Proposition B.3.

Proof B.4.

Lemma B.5.

Proof B.6.