Multilinear Time Invariant System Theory

TL;DR

This paper introduces a novel multilinear time invariant (MLTI) system framework using tensor algebra, extending classical LTI concepts to tensor-based dynamic systems for applications in biological and engineering contexts.

Contribution

It proposes a new MLTI system model based on Einstein product and paired tensors, extending classical system theory concepts to tensor representations.

Findings

Extended stability, reachability, and observability notions to MLTI systems.

Utilized recent tensor algebra advances for system analysis.

Provided a framework for tensor-based dynamic system modeling.

Abstract

In biological and engineering systems, structure, function and dynamics are highly coupled. Such interactions can be naturally and compactly captured via tensor based state space dynamic representations. However, such representations are not amenable to the standard system and controls framework which requires the state to be in the form of a vector. In order to address this limitation, recently a new class of multiway dynamical systems has been introduced in which the states, inputs and outputs are tensors. We propose a new form of multilinear time invariant (MLTI) systems based on the Einstein product and even-order paired tensors. We extend classical linear time invariant (LTI) system notions including stability, reachability and observability for the new MLTI system representation by leveraging recent advances in tensor algebra.

Equations

$$(\mathcal{X}\circ \mathcal{Y}{)}_{j_1{j}_2\dots j_N{i}_1{i}_2\dots i_M}={\mathcal{X}}_{j_1{j}_2\dots j_N}{\mathcal{Y}}_{i_1{i}_2\dots i_M}.$$

$$⟨\mathcal{X},\mathcal{Y}⟩=\sum _{\text{j}=\text{1}}^{\mathcal{J}}{\mathcal{X}}_{j_1{j}_2\dots j_N}{\mathcal{Y}}_{j_1{j}_2\dots j_N}$$

$$(\mathcal{X}{\times }_n\text{A}{)}_{j_1{j}_2\dots j_{n-1}i{j}_{n+1}\dots j_N}=\sum _{j_n=1}^{J_n}{\mathcal{X}}_{j_1{j}_2\dots j_n\dots j_N}{\text{A}}_{i{j}_n}.$$

$$\begin{array}{rl}& \mathcal{X}{\times }_1{\text{A}}_1{\times }_2\cdots {\times }_N{\text{A}}_N\\ & =\mathcal{X}\times \{{\text{A}}_1,{\text{A}}_2\dots ,{\text{A}}_N\}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N},\end{array}$$

$$ivec(\text{j},\mathcal{J})=j_1+\sum _{k=2}^N(j_k-1)\prod _{l=1}^{k-1}{J}_l.$$

$${\text{X}}_{\text{r}\times \text{c}}(j,i)={\mathcal{X}}_{j_1{j}_2\dots j_z{i}_1{i}_2\dots i_{N-z}}^{\mathbb{S}},$$

$$(\mathcal{A}\ast \mathcal{B}{)}_{j_1{i}_1\dots j_N{i}_N}=\sum _{\text{k}=\text{1}}^{\mathcal{K}}{\mathcal{A}}_{j_1{k}_1\dots j_N{k}_N}{\mathcal{B}}_{k_1{i}_1\dots k_N{i}_N}.$$

$${\mathcal{A}}_{j_1{i}_1\dots j_N{i}_N}\stackrel{\phi }{\to }{\text{A}}_{ivec(\text{j},\mathcal{J})ivec(\text{i},\mathcal{I})},$$

$$h(\mathcal{X})={\mathcal{X}}^{\top }\ast \mathcal{A}\ast \mathcal{X}>0$$

$$\begin{array}{rl}N(\mathcal{A})& =\{\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}:\mathcal{A}\ast \mathcal{X}=\mathcal{O}\},\\ R(\mathcal{A})& =\{\mathcal{A}\ast \mathcal{X}:\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}\},\end{array}$$

$$\displaystyle\text{rank}_{U}(\mathcal{A})=\text{rank}\big{(}\varphi(\mathcal{A})\big{)}.$$

$$({\Vert \begin{array}{cc}\mathcal{A}& \mathcal{B}\end{array}\Vert }_n{)}_{j_1{l}_1\dots j_N{l}_N}=\{\begin{array}{l}{\mathcal{A}}_{j_1{l}_1\dots j_N{l}_N},\text{j}\in [\mathcal{J}],\text{l}\in [\mathcal{I}]\\ {\mathcal{B}}_{j_1{l}_1\dots j_N{l}_N},\text{j}\in [\mathcal{J}],\text{l}\in [\mathcal{L}],\end{array}$$

$$\phi ({\Vert \begin{array}{cc}\mathcal{A}& \mathcal{B}\end{array}\Vert }_n)=\left[\begin{array}{cc}\phi (\mathcal{A})& \phi (\mathcal{B})\end{array}\right]\text{P},$$

$${\text{rank}}_U({\Vert \begin{array}{cc}\mathcal{A}& \mathcal{B}\end{array}\Vert }_n)=\text{rank}(\left[\begin{array}{cc}\phi (\mathcal{A})& \phi (\mathcal{B})\end{array}\right]).$$

$$\mathcal{A}\ast \mathcal{X}=\lambda \mathcal{X},$$

$$\mathcal{A}=\mathcal{V}\ast \mathcal{D}\ast {\mathcal{V}}^{-1},$$

$$\{\begin{array}{ll}{\mathcal{X}}_{t+1}& ={\mathcal{X}}_t\times \{{\text{A}}_1,\dots ,{\text{A}}_N\}+{\mathcal{U}}_t\times \{{\text{B}}_1,\dots ,{\text{B}}_N\}\\ {\mathcal{Y}}_t& ={\mathcal{X}}_t\times \{{\text{C}}_1,\dots ,{\text{C}}_N\}\end{array},$$

$$\{\begin{array}{ll}{\mathcal{X}}_{t+1}& =\mathcal{A}\ast {\mathcal{X}}_t+\mathcal{B}\ast {\mathcal{U}}_t\\ {\mathcal{Y}}_t& =\mathcal{C}\ast {\mathcal{X}}_t\end{array},$$

$${\mathcal{Y}}_{i_1{i}_2\dots i_N}$$

$$=\sum _{\text{j}=\text{1}}^{\mathcal{J}}{\mathcal{A}}_{i_1{j}_1\dots i_N{j}_N}{\mathcal{X}}_{j_1{j}_2\dots j_N},$$

$${\mathcal{X}}_{t+1}=\mathcal{A}\ast {\mathcal{X}}_t,$$

$${\mathcal{X}}_k={\mathcal{A}}^{\ast k}\ast {\mathcal{X}}_0+\sum _{j=0}^{k-1}{\mathcal{A}}^{\ast k-j-1}\ast \mathcal{B}\ast {\mathcal{U}}_j.$$

$${\mathcal{A}}^{\ast k}=\sum _{\text{j}=\text{1}}^{\mathcal{J}}{\lambda }_{j_1{j}_1\dots j_N{j}_N}^k{\mathcal{W}}_{j_1{j}_1\dots j_N{j}_N},$$

$${\mathcal{W}}_c(t_0,t_1)=\sum _{t=t_0}^{t_1-1}{\mathcal{A}}^{\ast t_1-t-1}\ast \mathcal{B}\ast {\mathcal{B}}^{\top }\ast ({\mathcal{A}}^{\top }{)}^{\ast t_1-t-1},$$

$${\mathcal{X}}_{t_1}$$

$$={\mathcal{A}}^{\ast t_1}\ast {\mathcal{X}}_0+{\mathcal{W}}_c(t_0,t_1)\ast {\mathcal{W}}_c^{-1}(t_0,t_1)\ast \mathcal{V}$$

$$={\mathcal{A}}^{\ast t_1}\ast {\mathcal{X}}_0+\mathcal{V}.$$

$${\mathcal{X}}_a+{\mathcal{A}}^{\ast t_1}\ast {\mathcal{X}}_0={\mathcal{A}}^{\ast t_1}\ast {\mathcal{X}}_0+\sum _{j=t_0}^{t_1-1}{\mathcal{A}}^{\ast t_1-j-1}\ast \mathcal{B}\ast {\mathcal{U}}_j.$$

$${\mathcal{X}}_a^{\top }\ast {\mathcal{X}}_a=\sum _{j=t_0}^{t_1-1}{\mathcal{X}}_a^{\top }\ast {\mathcal{A}}^{\ast t_1-j-1}\ast \mathcal{B}\ast {\mathcal{U}}_j=0,$$

$${\mathcal{W}}_c-\mathcal{A}\ast {\mathcal{W}}_c\ast {\mathcal{A}}^{\top }=\mathcal{B}\ast {\mathcal{B}}^{\top }.$$

$$\mathcal{R}=\Vert \begin{array}{cccc}\mathcal{B}& \mathcal{A}\ast \mathcal{B}& \dots & {\mathcal{A}}^{\ast |\mathcal{J}|-1}\ast \mathcal{B}\end{array}\Vert$$

$$\phi (\mathcal{R})=\left[\begin{array}{cccc}\text{B}& \text{A}\text{B}& \dots & {\text{A}}^{|\mathcal{J}|-1}\text{B}\end{array}\right]\text{P},$$

Full text

1 Introduction

In many complex systems, such as those arising in biology, capturing the interplay of function, structure and dynamics is critical in order to characterize the underlying mechanisms [1]. The human genome is a beautiful example of a multiway dynamical system [2]. The organization of the interphase nucleus reflects a dynamical interaction between 3D genome structure, function, and its relationship to phenotype, a concept known as the 4D Nucleome (4DN) [3]. 4DN research requires a comprehensive view of genome-wide structure, gene expression, the proteome, and phenotype which fits naturally with a tensorial representation [1]. The mathematical foundation of tensor based representation and analysis could play a critical role in the study of the human genome as well as social networks, cognitive science, signal processing and machine learning.

The notion of multilinear dynamical system or multilinear time invariant (MLTI) system was first introduced by Rogers et al. [4] for modeling of tensor time series, and Surana et al. [1] built the model by using tensor Tucker products to capture the evolution of the multilinear dynamics. Compared to classical linear time invariant (LTI) based approaches which fit vector or matrix models to tensor time series, MLTI representation provides a more natural, compact and accurate representation of tensorial data with fewer model parameters. By using tensor unfolding, an operation that transforms a tensor into a matrix, Rogers et al. [4] and Surana et al. [1] developed methods for model identification/reduction from tensor time series data, and demonstrated benefits of the MTLI representation compared to the classical LTI approach. However, this representation of MLTI systems is limited by the fact that it assumes the multilinear operators are formed from the Tucker products of matrices, and thus precludes more general tensorial representation. Moreover, while tensor unfolding enables one to transform MLTI into a classical LTI representation for computational purposes (i.e. enables use of matrix algebra), one loses the inherent tensor algebraic structure which otherwise could be exploited to develop system theoretic concepts.

In this paper, we propose a more generalized MLTI system model using the Einstein product, from which the Tucker product based MTLI representation can be obtained as a special case. The Einstein product is a tensor contraction operation used quite often in tensor calculus and has profound applications in the study of continuum mechanics and the field of relativity theory [5, 6]. The proposed generalized MLTI system model takes a very similar form to the classical LTI system model, and is thus more naturally suited to develop system theoretic concepts. Moreover, the space of even-order tensors equipped with the Einstein product has many desirable properties. Brazell et al. [7] in 2013 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. group of invertible matrices.

Building on the results of Brazell et al. [7] and the notion of block tensors [8], we propose new tensor notions for positive definiteness, unfolding rank and a new way of concatenation of tensors to create block tensors. Using these tensor constructs, we develop tensor algebraic conditions for stability, reachability and observability for the generalized MLTI systems. Interestingly, these new conditions look analogous to the classical conditions for stability, reachability and observability for LTI systems, and reduce to them in special cases. To the best of the authors’ knowledge, MTLI system representation using the Einstein product and formulation of these tensor algebraic system theoretic conditions has never been reported in the literature. Due to space limitations, we only provide proofs for selected results. A more comprehensive publication is under preparation.

The paper is organized in five sections. We start with basics of tensor algebra followed with tensor groups, block tensors and tensor eigenvalue decompositions in Section 2. A new general representation of MLTI systems is introduced in Section 3, and generalization of stability, reachability and observability conditions for the MLTI systems is also discussed. A simple single input single output MLTI system example is given in Section 4, and we conclude with directions for future research in Section 5.

2 Tensor Algebra

We use concepts and notations for tensor algebra from the comprehensive works of Kolda et al. [9, 10] and Ragnarsson et al. [8]. A tensor is a multidimensional array. The order of a tensor is the number of its dimensions. An $N$-th order tensor usually is denoted by $\mathcal{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$. The sets of indexed indices and size of $\mathcal{X}$ are denoted by $\textbf{j}=\{j_{1},j_{2},\dots,j_{N}\}$ and $\mathscr{J}=\{J_{1},J_{2},\dots,J_{N}\}$, respectively. $|\mathscr{J}|$ represents the product of all elements in $\mathscr{J}$, and $\textbf{j}\in[\mathscr{J}]$ can be interpreted as $j_{n}=1,2,\dots,J_{n}$ for $n=1,2,\dots,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.

By extending the notion of vector outer product, the outer product of two tensors $\mathcal{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ and $\mathcal{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 $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ is defined as

[TABLE]

leading to the Frobenius norm $\|\mathcal{X}\|^{2}=\langle\mathcal{X},\mathcal{X}\rangle$. The notation $\sum_{\textbf{j}=\textbf{1}}^{\mathscr{J}}$ can be read as an abbreviation of $N$ summations over all indices $j_{n}=1,2,\dots,J_{n}$ for $n=1,2,\dots,N$.

The matrix tensor multiplication $\mathcal{X}\times_{n}\textbf{A}$ along mode $n$ for a matrix $\textbf{A}\in\mathbb{R}^{I\times J_{n}}$ is defined by

[TABLE]

This product can be generalized to what is known as the Tucker product,

[TABLE]

where, $\textbf{A}_{n}\in\mathbb{R}^{I_{n}\times J_{n}}$. If an $N$-th order tensor $\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{N}}$ can be expressed as $\mathcal{Y}=\mathcal{X}\times\{\textbf{A}_{1},\textbf{A}_{2}\dots,\textbf{A}_{N}\}$, the decomposition is referred to as the Tucker decomposition. When the core tensor $\mathcal{X}$ possesses the “higher-order diagonal” property and the factor matrices $\textbf{A}_{n},n=1,\cdots,N$ are unitary, it is also called the Higher-Order Singular Value Decomposition (HOSVD), a multilinear generalization of the matrix Singular Value Decomposition (SVD) [11].

Tensor unfolding is considered as a critical operation in tensor computations [8, 9, 10]. In order to unfold a tensor $\mathcal{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ into a matrix, we use an index mapping function $ivec(\cdot,\mathscr{J}):\underbrace{\mathbb{Z}^{+}\times\mathbb{Z}^{+}\times\dots\times\mathbb{Z}^{+}}_{N}\rightarrow\mathbb{Z}^{+}$ defined by Ragnarsson et al. [8], which is given as

[TABLE]

Suppose that $\textbf{r}=\mathbb{S}(1:z)$ and $\textbf{c}=\mathbb{S}(z+1:N)$, where $z$ is an integer such that $1\leq z<N$, and $\mathbb{S}$ is a vector from the set of all permutations of $1$ to $N$. Define $\mathscr{J}(\textbf{r})=\{J_{\mathbb{S}(1)},J_{\mathbb{S}(2)},\dots,J_{\mathbb{S}(z)}\}$ and $\mathscr{J}(\textbf{c})=\{J_{\mathbb{S}(z+1)},J_{\mathbb{S}(z+2)},\dots,J_{\mathbb{S}(N)}\}$. Then the $\textbf{r}\times\textbf{c}$ unfolding matrix of $\mathcal{X}$ denoted by $\textbf{X}_{\textbf{r}\times\textbf{c}}$ is given by

[TABLE]

where, $j=ivec(\textbf{j},\mathscr{J}(\textbf{r}))$, $i=ivec(\textbf{i},\mathscr{J}(\textbf{c}))$ and $\mathcal{X}^{\mathbb{S}}$ is the $\mathbb{S}$-transpose of $\mathcal{X}$ (see (2.9) in [8]). In particular, when $z=1$ and $\mathbb{S}=\{n,1:n-1,n+1:N\}$, the tensor unfolding is called the n-mode matricization.

2.1 Einstein Product and Isomorphism

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. Similarly, these notions proved to be useful in numerical solutions of master equations associated with Markov processes on extremely large state spaces [12]. Even-order paired tensors were originally proposed by Huang and Qi [13] in the context of elasticity tensors in solid mechanics. It turns out that compared to even-order non-paired tensors, the even-order paired tensors are easier for bookkeeping and can be conveniently manipulated using tensor algebra for MLTI systems, see Section 3.

Definition 2.1

Given an even-order tensor $\mathcal{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$, if its indices can be divided into $N$ adjacent blocks $\{j_{1}i_{1}\},\dots,\{j_{N}i_{N}\}$, then $\mathcal{A}$ is called an even-order paired tensor.

Definition 2.2

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

[TABLE]

If $\mathcal{B}\in\mathbb{R}^{K_{1}\times K_{2}\times\dots\times K_{N}}$, the Einstein product is still valid by treating $\mathscr{I}=\textbf{1}$ in (2.7). Note that the notion of Einstein product is not restricted to even-order paired tensors and can be defined more generally.

The above Einstein product computes the summation of two even-order paired tensors over alternating indices. Brazell et al. [7] investigated properties for even-order non-paired tensors $\mathcal{A}\in\mathbb{R}^{J_{1}\times\dots\times J_{N}\times I_{1}\times\dots\times I_{N}}$ under the Einstein product through construction of an isomorphism to GL($\mathbb{R}$) (general linear group). The existence of the isomorphism enables one to generalize several matrix concepts, such as orthogonality, invertibility and eigenvalue decomposition to the tensor case [7, 14, 15, 16, 17]. We establish an analogous isomorphism for even-order paired tensors by a permutation of indices.

Definition 2.3

Define the following transformation $\varphi$: $\mathbb{T}_{J_{1}I_{1}\dots J_{N}I_{N}}(\mathbb{R})\rightarrow\mathbb{M}_{|\mathscr{J}|,|\mathscr{I}|}(\mathbb{R})$ with $\varphi(\mathcal{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 $J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}$ even-order paired tensors and $\mathbb{M}_{|\mathscr{J}|,|\mathscr{I}|}(\mathbb{R})$ is the set of all $|\mathscr{J}|\times|\mathscr{I}|$ matrices.

The transformation $\varphi$ can be viewed as a tensor unfolding discussed in (2.6) with $z=N$ and $\mathbb{S}=\{1,3,\dots,2N-1,2,4,\dots,2N\}$. Brazell et al. prove that $\varphi$ is an isomorphism for fourth-order non-paired tensors, and we extend the key results (Corollary 3.3 in [7]) for even-order paired tensors of any order.

Corollary 2.1

Suppose that $J_{n}=I_{n}$ for all $n$ and $\mathbb{M}_{|\mathscr{J}|,|\mathscr{I}|}(\mathbb{R})=\text{GL}(|\mathscr{J}|,\mathbb{R})$. $\mathbb{T}_{J_{1}J_{1}\dots J_{N}J_{N}}(\mathbb{R})$ is a group equipped with the Einstein product, and $\varphi$ is a group isomorphism. Moreover, $\mathbb{T}_{J_{1}J_{1}\dots J_{N}J_{N}}(\mathbb{R})$ also forms a tensor ring under addition and the Einstein product, and $\varphi$ is a ring isomorphism.

For an even-order paired tensor $\mathcal{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$, $\mathcal{T}\in\mathbb{R}^{I_{1}\times J_{1}\times\dots\times I_{N}\times J_{N}}$ is called the U-transpose of $\mathcal{A}$ if $\mathcal{T}_{i_{1}j_{1}\dots i_{N}j_{N}}=\mathcal{A}_{j_{1}i_{1}\dots j_{N}i_{N}}$ and is denoted by $\mathcal{A}^{\top}$. We refer to an even-order paired tensor that is identical to its U-transpose as weakly symmetric. An even-order “square” tensor $\mathcal{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 $\mathcal{D}_{j_{1}j_{1}\dots j_{N}j_{N}}$. In particular, if all the diagonal entires $\mathcal{D}_{j_{1}j_{1}\dots j_{N}j_{N}}=1$, then $\mathcal{D}$ is the U-identity tensor, denoted by $\mathcal{I}$. An even-order square tensor $\mathcal{U}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is U-orthogonal if $\mathcal{U}*\mathcal{U}^{\top}=\mathcal{U}^{\top}*\mathcal{U}=\mathcal{I}$. Furthermore, for an even-order square tensor $\mathcal{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$, if there exists a tensor $\mathcal{X}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ such that $\mathcal{A}*\mathcal{X}=\mathcal{X}*\mathcal{A}=\mathcal{I}$, then $\mathcal{X}$ is called the U-inverse of $\mathcal{A}$, denoted by $\mathcal{A}^{-1}$. “U” stands for “unfolding” in all the definitions. Besides the properties discussed above, we can define U-positive definiteness for even-order square tensors, which is similarly discussed in [12].

Definition 2.4

An even-order square tensor $\mathcal{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is U-positive definite if its corresponding homogeneous polynomial

[TABLE]

for all $\mathcal{X}\neq\mathcal{O}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$, where $\mathcal{O}$ denotes the zero tensor.

It is straightforward to show that an even-order square tensor $\mathcal{A}$ is U-positive definite if and only if $\varphi(\mathcal{A})$ is positive definite. Moreover, U-positive definiteness implies U-invertibility of even-order square tensors from the isomorphism property.

The notions of linear dependence and independence for tensor spaces are defined in a similar way to those for vector spaces. Let $\mathbb{T}_{J_{1}J_{2}\dots J_{N}}(\mathbb{R})$ be a set of all tensors $\mathcal{X}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$. A basis $\mathscr{B}$ of the tensor space $\mathbb{T}_{J_{1}J_{2}\dots J_{N}}(\mathbb{R})$ is a linearly independent subset of $\mathbb{T}_{J_{1}J_{2}\dots J_{N}}(\mathbb{R})$ that spans $\mathbb{T}_{J_{1}J_{2}\dots J_{N}}(\mathbb{R})$. Additionally, \text{dim}\big{(}\mathbb{T}_{J_{1}J_{2}\dots J_{N}}(\mathbb{R})\big{)} is equal to the cardinality of the basis $\mathscr{B}$. Ji et al. [14] propose the null space $N(\mathcal{A})$ and the range $R(\mathcal{A})$ of an even-order non-paired tensor $\mathcal{A}$, and establish some elementary properties. Analogous definitions for even-order paired tensors are as follows:

Definition 2.5

Define the null space and range of an even-order paired tensor $\mathcal{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ to be

[TABLE]

respectively. Moreover, define nullity${}_{U}(\mathcal{A}$)=dim($N(\mathcal{A})$) and rank${}_{U}(\mathcal{A}$)=dim($R(\mathcal{A})$).

We find that above notion of rank${}_{U}(\mathcal{A}$) is equivalent to the unfolding rank defined by Liang et al. [15].

Proposition 2.1

Let $\mathcal{A}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be an even-order paired tensor. Then

[TABLE]

The proof is based on the definitions of the unfolding rank and the bases in tensor spaces as discussed above. Hence, we also refer to the rank${}_{U}(\mathcal{A}$) as the unfolding rank in this paper. In addition, Liang et al. [15] propose the unfolding determinant, $\text{det}_{U}$, for even-order non-paired tensors. Similar definition/results can be extended for even-order paired tensors.

2.2 Block Tensors

Analogously to block matrices, one can define the notion of block tensors. The block tensors introduced by Sun et al. [16] for even-order non-paired tensors have the limitation of introducing too many zeros into the tensor and increasing its size. For tensors of the same size, we explore a new block tensor construction that does not introduce any wasteful zeros, and thus could offer computational advantage.

Definition 2.6

Let $\mathcal{A},\mathcal{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be two even-order paired tensors of the same size. Then the $n$-mode row block tensor is defined to be $\begin{Vmatrix}\mathcal{A}&\mathcal{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}}$ such that

[TABLE]

where, $\mathscr{L}=\mathscr{I}$ except $l_{n}=I_{n}+1,I_{n}+2,\dots,2I_{n}$.

The $n$-mode column block tensor $\begin{Vmatrix}\mathcal{A}&\mathcal{B}\end{Vmatrix}_{n}^{\top}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times 2J_{n}\times I_{n}\times\dots\times J_{N}\times I_{N}}$ can be defined in a similar manner. However, the blocks of even-order paired tensors usually do not map to contiguous blocks in their unfolding [8], and sometimes that may increase the complexity in computations. Ragnarsson et al. [8] show that there exists a row permutation matrix Q and a column permutation matrix P such that the blocks of a tensor $\mathcal{A}$ can be mapped to contiguous blocks in the unfolding $\textbf{Q}\textbf{A}_{\textbf{r}\times\textbf{c}}\textbf{P}$. The following proposition shows that for $n$-mode row block tensors, only column permutations are required.

Proposition 2.2

Let $\mathcal{A},\mathcal{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be two even-order paired tensors. Then

[TABLE]

where, P is a permutation matrix. In particular, when $\mathscr{I}=\textbf{1}$ or $n=N$, P is the identity matrix.

The proof follows immediately from the definition of index mapping function $ivec(\textbf{i},\mathscr{I})$ and $n$-mode row block tensors. Based on Ragnarsson et al.’s results, it follows that the blocks of $n$-mode column block tensors map to contiguous blocks in its unfolding up to some row permutations. Moreover, Proposition 2.2 helps us to establish the relation between unfolding rank and matrix rank for block tensors and their unfolding.

Corollary 2.2

Let $\mathcal{A},\mathcal{B}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$ be two even-order paired tensors. Then

[TABLE]

Lastly, we generalize the $n$-mode block tensors for multiple blocks. Given $K$ even-order paired tensors $\mathcal{X}_{n}\in\mathbb{R}^{J_{1}\times I_{1}\times\dots\times J_{N}\times I_{N}}$, one can apply Definition 2.6 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. We define a more general concatenation approach as follows:

Definition 2.7

Given $K$ even-order paired tensors $\mathcal{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 $\mathcal{Y}$ can be constructed in the following way:

• •

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

• •

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

• •

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

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

The generalized mode column block tensors with multiple blocks can be constructed in a similar manner. Both Proposition 2.2 and Corollary 2.2 hold for the general mode block tensor construction defined above. In particular, when $\mathscr{I}=\textbf{1}$, 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, e.g. see next subsection.

2.3 Tensor Eigenvalue Decomposition

Eigenvalue problems for higher-order tensors were first explored by Qi [18] and Lim [19] independently in 2005. Brazell et al. [7] formulated a new tensor eigenvalue problem through the isomorphism $\varphi$ for fourth-order non-paired tensors, and Cui et al. [17] extended the result to general even-order tensors. Note that Tensor Eigenvalue Decomposition (TEVD) derived from the unfolding transformation is distinct from other notions of tensor decompositions, such as the Candecomp/Parafac decomposition (CP decomposition) or the Tucker decomposition. We exploit our new notion of mode block tensors to express the full tensor eigenvalue decomposition of an even-order paired tensor.

Definition 2.8

Let $\mathcal{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ be an even-order square tensor. If $\mathcal{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 $\mathcal{X}$ and $\lambda$ satisfy

[TABLE]

then we call $\lambda$ and $\mathcal{X}$ as the U-eigenvalue and U-eigentensor of $\mathcal{A}$, respectively. Moreover, the tensor eigenvalue decomposition is given by

[TABLE]

where, $\mathcal{D}$ is an U-diagonal tensor with U-eigenvalues on its diagonal, and $\mathcal{V}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is a mode row block tensors consisting of all the U-eigentensors, i.e. $\mathcal{V}=\begin{Vmatrix}\mathcal{X}_{1}&\mathcal{X}_{2}&\dots&\mathcal{X}_{|\mathscr{J}|}\end{Vmatrix}$. Here we have chosen $K_{1}=J_{1},\cdots,K_{n}=J_{n}$ in applying the mode row block tensor operation which enables to express the TEVD in the form (2.16) analogous to the matrix case.

The algebraic and geometric multiplicity of tensor eigenvalues can be defined as for matrices. The following provides a generalization of the Caley-Hamilton theorem for the tensor case.

Lemma 2.1

If $\mathcal{A}\in\mathbb{R}^{J_{1}\times J_{1}\times\dots\times J_{N}\times J_{N}}$ is an even-order square tensor, then $\mathcal{A}$ satisfies its own characteristic polynomial $p(\lambda)=\text{det}_{U}(\lambda\mathcal{I}-\mathcal{A})$, i.e. $p(\mathcal{A})=\mathcal{O}$.

• Proof.

The proof follows immediately by applying $\varphi$.

3 MLTI Systems Theory

In order to describe the evolution of tensor time series, the authors in [1, 4] introduced a MLTI system involving the Tucker product as follows,

[TABLE]

where, $\mathcal{X}_{t}\in\mathbb{R}^{J_{1}\times J_{2}\times\dots\times J_{N}}$ is the latent state space tensor, $\mathcal{Y}_{t}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{N}}$ is the output tensor, and $\mathcal{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. In particular, using the Kronecker product, one can transform the system (3.17) into a standard LTI system [4] and then apply the standard LTI systems concepts for analysis. However, this representation is limited by several factors. Firstly, the set of multilinear operators resulting from the component matrices $\textbf{A}_{n}$, $\textbf{B}_{n}$ and $\textbf{C}_{n}$ consists of a special case and does not capture the more general multilinear evolution of tensor dynamics (see system (3.18)). Secondly, once transformed into an LTI system via the Kronecker product, there is no unique way to recover the original tensor based representation. Thus, one loses the inherent tensor algebraic structure which otherwise could be exploited to develop system theoretic concepts such as reachability and observability Gramians more naturally.

We find that (3.17) can be replaced by a more general representation using the notion of even-order paired tensors and the Einstein product.

Definition 3.1

A more general representation of MLTI system is given by

[TABLE]

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

Proposition 3.1

The governing equations (3.18) can be obtained from (3.17) by setting $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{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.

Consider a general Tucker decomposition $\mathcal{Y}=\mathcal{X}\times\{\textbf{A}_{1},\textbf{A}_{2},\dots,\textbf{A}_{N}\}$ and rewrite it elementwise, i.e.

[TABLE]

for $\mathcal{A}_{i_{1}j_{1}\dots i_{N}j_{N}}=(\textbf{A}_{1})_{i_{1}j_{1}}\dots(\textbf{A}_{N})_{i_{N}j_{N}}$. By the definitions of outer product and the Einstein product, it follows that $\mathcal{Y}=\mathcal{A}*\mathcal{X}$ where $\mathcal{A}=\textbf{A}_{1}\circ\textbf{A}_{2}\circ\dots\circ\textbf{A}_{N}$. Hence, the result follows immediately.

The Einstein product representation (3.18) of MLTI systems is indeed the generalization of (3.17), and overcomes most of the limitations of Tucker product based MTLI representation discussed above. More importantly, it takes a form similar to the standard LTI system model, and so the representation is more natural for developing MLTI systems theory which we discuss next.

3.1 Solution of MLTI System

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

Proposition 3.2

For the unforced MLTI system

[TABLE]

the solution for $\mathcal{X}$ at time $k$ given initial condition $\mathcal{X}_{0}$ is $\mathcal{X}_{k}=\mathcal{A}^{*k}*\mathcal{X}_{0}$ where $\mathcal{A}^{*k}=\underbrace{\mathcal{A}*\mathcal{A}*\dots*\mathcal{A}}_{k}$.

The proof is straightforward using the notion of even-order paired tensors and the Einstein product. If the even-order paired tensor is of the form $\mathcal{A}=\textbf{A}_{1}\circ\textbf{A}_{2}\circ\dots\circ\textbf{A}_{N}$, the $k$-th power Einstein product of $\mathcal{A}$ can be computed by $\mathcal{A}^{*k}=\textbf{A}_{1}^{k}\circ\textbf{A}_{2}^{k}\circ\dots\circ\textbf{A}_{N}^{k}$. Applying Proposition 3.2, we can write down the explicit solution of (3.18) which takes an analogous form to the LTI system,

[TABLE]

3.2 Stability.

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

Proposition 3.3

For an unforced MLTI system (3.19), the equilibrium point $\mathcal{X}=\mathcal{O}$ is:

• •

stable if and only if the magnitudes of all the U-eigenvalues of $\mathcal{A}$ are less than or equal to 1; for those equal to 1, its algebraic and geometry multiplicity must be equal;

• •

asymptotically (or exponentially) stable if the magnitudes of all the U-eigenvalues are less than 1;

• •

unstable if the magnitudes of some of the U-eigenvalues are greater than 1.

• Proof.

We only focus on the case when $\mathcal{A}$ has a full set of U-eigentensors, i.e. $\mathcal{A}=\mathcal{V}*\mathcal{D}*\mathcal{V}^{-1}$ in (2.16). It follows from the solution of system (3.19) that

[TABLE]

where, $\lambda_{j_{1}j_{1}\dots j_{N}j_{N}}$ are the U-eigenvalues of $\mathcal{A}$, and $\mathcal{W}_{j_{1}j_{1}\dots j_{N}j_{N}}$ are some even-order paired tensors. Then the results follow immediately.

3.3 Reachability

Here and in the following subsection, we introduce the definitions of reachability and observability for MLTI systems which are similar to analogous concepts for the LTI systems [20, 21, 22]. We then establish sufficient and necessary conditions for reachability and observability for the MLTI systems.

Definition 3.2

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

Theorem 3.1

The pair $(\mathcal{A},\mathcal{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.

Suppose $\mathcal{W}_{c}(t_{0},t_{1})$ is U-positive definite, and let $\mathcal{X}_{0}$ be the initial state and $\mathcal{X}_{1}$ be the desired final state. Choose $\mathcal{U}_{t}=\mathcal{B}^{\top}*(\mathcal{A}^{\top})^{*t_{1}-t-1}*\mathcal{W}_{c}^{-1}(t_{0},t_{1})*\mathcal{V}$ for some constant tensor $\mathcal{V}$. It follows from the explicit solution of system (3.18) that

[TABLE]

Take $\mathcal{V}=-\mathcal{A}^{*t_{1}}*\mathcal{X}_{0}+\mathcal{X}_{1}$, we have $\mathcal{X}_{t_{1}}=\mathcal{X}_{1}$.

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

[TABLE]

Multiplying from the left by $\mathcal{X}_{a}^{\top}$ yields

[TABLE]

which implies that $\mathcal{X}_{a}=\mathcal{O}$, a contradiction.

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

[TABLE]

If the pair $(\mathcal{A},\mathcal{B})$ is reachable and all the U-eigenvalues of $\mathcal{A}$ have magnitude less than 1, one can show that there exists a unique weakly symmetric U-positive definite solution $\mathcal{W}_{c}$. 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 3.4

The pair $(\mathcal{A},\mathcal{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})=|\mathscr{J}|$.

• Proof.

Using the unfolding transformation $\varphi$ one can express,

[TABLE]

where, $\textbf{A}=\varphi(\mathcal{A})$, $\textbf{B}=\varphi(\textbf{B})$ and P is some permutation matrix. Then the result follows immediately.

Alternatively, the proof can be developed directly in the tensor setting by applying Lemma 2.1 like the approach used in Theorem 3.1. When $N=1$, Corollary 3.4 simplifies to the famous Kalman rank condition for reachability of LTI systems.

3.4 Observability

Definition 3.3

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

Theorem 3.2

The pair $(\mathcal{A},\mathcal{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.

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

[TABLE]

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

Proposition 3.5

The pair $(\mathcal{A},\mathcal{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})=|\mathscr{J}|$.

4 Numerical Example

To illustrate MLTI systems theory, we consider a simple single input single output system that is given by (3.17) with

[TABLE]

and the states $\mathcal{X}_{t}\in\mathbb{R}^{3\times 2}$ are second-order tensors, i.e. matrices. The U-eigenvalues of $\mathcal{A}=\textbf{A}_{1}\circ\textbf{A}_{2}$ are $\pm 0.9207$, $-0.1775\pm 0.2128i$ and $0.1775\pm 0.2128i$. Hence, the MLTI system is asymptotically stable. In addition, the reachability and observability tensors are given by

[TABLE]

and

[TABLE]

Furthermore, $\text{rank}_{U}(\mathscr{R})=6$ and $\text{rank}_{U}(\mathscr{O})=6$, and the system therefore is both reachable and observable. For all the computations in this example, we used the unfolding transform $\varphi$ which enabled us to use standard matrix algebra.

Remark: Note that unfolding transform allows one to transform tensor algebra problems to standard matrix algebra problems. However, it may not be the most memory and numerically efficient approach. In fact, computing tensor algebraic notions without unfolding is an active area of research [7, 12, 15, 24], and we are currently exploring methods based on that for computations associated with the MLTI systems.

5 Conclusion

In this paper, we generalized the MLTI system representation using even-order paired tensors and the Einstein product. This new representation also facilitated the generalization of notions of stability, reachability and observability from classical multivariate control theory to that for MLTI systems. In particular, the unfolding isomorphism played a key role in establishing criterion for stability, reachability and observability for MLTI systems.

In future work, it should be worthwhile to develop an associated theoretical and computational framework for data driven model identification/reduction, observer and feedback control design, and to apply these techniques to real world engineering systems and machine learning. 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 [25, 26]. Such applications would also ideally be analyzed using nonlinearity and stochasticity in the multiway dynamical system representation and analysis framework. This is an important direction for future research.

Acknowledgments: We thank Dr. Frederick Leve at the Air Force Office of Scientific Research (AFOSR) for support and encouragement. This work is supported in part under AFOSR Award No: FA9550-18-1-0028, NSF grant DMS 1613819 and the Lifelong Learning Machines program from DARPA/MTO.