{'fraction_non_alphanumeric': 0.05136179775280899, 'fraction_numerical': 0.013519101123595506, 'mean_word_length': 4.092090809227389, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 50, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-Driven Model Reduction for Multilinear Control Systems via Tensor Trains

In this paper, we explore the role of tensor algebra in balanced truncation (BT) based model reduction/identification for high-dimensional multilinear/linear time invariant systems. In particular, we employ tensor train decomposition (TTD), which provides a good compromise between numerical stability and level of compression, and has an associated algebra that facilitates computations. Using TTD, we propose a new BT approach which we refer to as higher-order balanced truncation, and consider different data-driven variations including higher-order empirical gramians, higher-order balanced proper orthogonal decomposition and a higher-order eigensystem realization algorithm. We perform computational and memory complexity analysis for these different flavors of TTD based BT methods, and compare with the corresponding standard BT methods in order to develop insights into where the proposed framework may be beneficial. We provide numerical results on simulated and experimental datasets showing the efficacy of the proposed framework.

1. Introduction. The goal of this paper is to explore the role of tensor algebra in dataguided model reduction/identification for high-dimensional multilinear/linear time invariant (MLTI/LTI) input-output systems. 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 computations and signal processing [7,9,11,19,27,28,29,62]. Tensor representation preserves multidimensional patterns, capturing higher-order interactions and couplings within multiway data, instead of the standard pairwise "flattened" view inherent in two-way matrix based analysis. Tensor decomposition techniques such as CANDECOMP/PARAFAC decomposition (CPD) [1,27,28,49], higherorder singular value decomposition (HOSVD) [12] and tensor train decomposition (TTD) [43,46] help reveal such hidden patterns/redundancies to obtain a compact representation, thereby reducing storage effort and enabling efficient computations.
Tensor algebra has been recently exploited in systems and control applications. The key idea is to tensorize the vector based dynamic system representation into an equivalent tensor representation, and to exploit tensor algebra. A large body of literature, referred to as tensor product models, has emerged which utilizes compact tensor based representations/computa-tions in the context of linear parameter varying models [3]. Other applications include efficient solution of Lyapunov equations [39], fault detection [38], Kalman filtering [4], accelerating simulation of nonlinear models where the vector field is a multilinear function of states [30] and modeling inverse dynamics [2], to name a few.
In many scientific and engineering applications, the system dynamics over space and time is often described in terms of partial differential equations (PDEs), e.g., Navier Stokes and heat equations in thermal/fluids, Euler equations in structural mechanics, Schrödinger equations in quantum mechanics, etc. The system state in such representation is a field (e.g., temperature/velocity field) whose discretization in space/time naturally results in tensors evolving with time. In order to apply the standard model reduction/identification framework, such as Balanced Truncation (BT) [54], Eigensystem Realization Algorithm (ERA) [24,36] or Dynamic Mode Decomposition (DMD) and its variants [50,58], tensors need to be vectorized. This may result in an extremely high-dimensional system representation in which the number of states/model parameters scales exponentially with the number of dimensions of the tensors involved, and thus pose a significant computational challenge. Alternatively, a new class of multilinear time invariant (MLTI) systems has been introduced [8,9,53,57] in which the states and outputs are preserved as tensors, and the system evolution is generated by the action of multilinear operators. By using tensor unfolding, an operation that transforms a tensor into a matrix, Rogers et al. [53] and Surana et al. [57] developed methods for model reduction/identification from tensor time series data, and demonstrated benefits such as a more compact and accurate representation compared to the classical vectorization based LTI approach. An application of such tensor based representation and identification for skeleton based human behavior recognition from videos demonstrated significant improvements in classification accuracy compared to standard LTI based approaches [13]. The TTD framework has been applied to accelerate DMD computations [26], and for computing numerical solutions of master equations associated with high-dimensional Markov processes [16].
Chen et al. [8,9] generalized the notion of MTLI systems to incorporate control inputs based on the Einstein product and even-order paired tensors. By leveraging recent advances in tensor algebra, Chen et al. also developed tensor generalizations of classical LTI system notions including stability, reachability and observability, and expressed them in different forms including concepts based on tensor unfolding and other more standard notions of tensor ranks/decompositions [9]. In addition, the authors in [9] demonstrated that the generalized CPD and TTD (GTTD) based model reduction framework can significantly reduce the number of MLTI system parameters. Continuing along that line of work, in this paper, we develop a novel tensor based BT computational framework for model reduction/identification in highdimensional MLTI/LTI systems. In particular, we choose to work with TTD as the underlying computational framework since it provides a good compromise between numerical stability and level of compression. The proposed framework is naturally suited for MLTI systems, but can also be applied to standard vector based LTI system representations by tensorizing them into a suitable tensor form [26,39,44]. The key contributions of this paper are as follows: 1. We develop a TTD based computational framework for different flavors of BT: Higher-Order Balanced Truncation (HOBT) which involves solutions of tensor Lyapunov equations to obtain reachability and observability gramians, and subsequent computations of the balancing transform; Higher-Order Empirical Gramians (HOEG) which constructs reachability and observability gramians from snapshots obtained via simulation of forward/adjoint MTLI system equations; Higher-Order Balanced Proper Orthogonal Decomposition (HOBPOD) which directly and more efficiently constructs the balancing transform from forward/adjoint snapshots; and Higher-Order Eigensystem Realization Algorithm (HOERA) which is equivalent to HOBPOD but only requires snapshots from the forward model, and so can also be applied to experimental data. 2. We employ TT-algebra for the memory efficient computations whereby basic tensor operations such as addition, scalar product, Einstein product, norms, solution to multilinear equations and tensor pseudoinverse, can be computed and maintained in the TT-format, without needing to go to the full tensor representation. We also present variations for computing block tensors and tensor singular value decomposition (TSVD) using TT-algebra, which are the key operations in HOBPOD/HOERA. 3. We provide computational and memory complexity analysis for different tensor based BT methods, and compare them with the corresponding standard BT methods. This analysis provides insights into when the tensor based representation can be beneficial, and also suggests a framework whereby one can mix tensor and matrix approaches to gain the best efficiency. 4. We demonstrate our framework in four numerical examples including a synthetic dataset, the 2D heat equation with control, a cancer cell image video dataset, and a room impulse responses dataset. A comparison with the standard BT is provided indicating the memory and computational time savings. The paper is organized into six sections. We start with the basics of tensor algebra followed by descriptions of even-order paired tensors, tensor decompositions and TT-algebra in section 2. We also establish results about the construction of block tensors and the computation of TSVD using TTD. In section 3, we present a generalization of the balanced truncation, balanced proper orthogonal decomposition and eigensystem realization algorithm framework for MLTI systems with detailed algorithms and complexity analysis. Four numerical examples are presented in section 4. Finally, we discuss some directions for future research in section 5 and conclude in section 6. For ease of reading, we provide a list of acronyms in Appendix B.

Preliminaries.
In this section, we briefly review tensor preliminaries. Comprehensive reviews can be found in [9,16,27,28,43]. An N -th order tensor usually is denoted by X ∈ R J 1 ×J 2 ×···×J N . The sets of indexed indices and size of X are denoted by j = {j 1 , j 2 , . . . , j N } and J = {J 1 , J 2 , . . . , J N }, respectively. |J | represents the product of all elements in J , and j ∈ [J ] can be interpreted as j n = 1, 2, . . . , J n for n = 1, 2, . . . , N . For more compact notation, we define R J 1 ×J 2 ×···×J N = R J and X j 1 j 2 ...j N = X j . The tensor inner product of two tensors of the same size is defined by X, Y = J j=1 X j Y j where the notation J j=1 is an abbreviation of the N summations over all indices j ∈ [J ]. The tensor Frobenius norm induced by the inner product is given as X 2 = X, X . The matrix tensor multiplication X × n A along mode n for a matrix A ∈ R I×Jn is defined by (X × n A) j 1 j 2 ...j n−1 ij n+1 ...j N = Jn jn=1 X j A ijn . This product can be generalized to what is known as the Tucker product, for A n ∈ R In×Jn , 2.1. Even-order paired tensors. The notion of even-order paired tensors was first proposed by Huang et al. [21] for elasticity tensors in the context of solid mechanics. For an arbitrary even-order tensor A ∈ R J 1 ×I 1 ×···×J N ×I N , if its indices can be divided into N adjacent blocks {j 1 i 1 }, . . . , {j N i N }, then A is called an even-order paired tensor [8]. Similarly, we define R J 1 ×I 1 ×···×J N ×I N = R J ⊗I and A j 1 i 1 ...j N i N = A j⊗i for simplicity. Building on the work by Brazell et al. [7], Chen et al. [8] proposed an unfolding transformation ψ from the even-order paired tensor space T J ⊗I (R) to the matrix space M |J ||I| (R) defined by and showed that ψ is a ring isomorphism for J = I under the element-wise addition and the Einstein product defined by for A ∈ R J ⊗K and B ∈ R K⊗I . Based on the transformation ψ, one can define matrix-like notions for tensor algebra including U-transpose, U-diagonal, U-identity, U-orthogonal, Uinverse, U-positive definiteness and U-eigenvalues (see [9] for details). Liang et al. [35] also define the unfolding rank of an even-order paired tensor A ∈ R J ⊗I by rank U (A) = rank ψ(A) . Chen et al. [8] introduce notions of block tensors for same-size even-order paired tensors. For A, B ∈ R J ⊗I , the n-mode row block tensor is defined to be A B n ∈ R J ⊗L such that where, L = {I 1 , I 2 , . . . , 2I n , . . . , I N } and l n = I n + 1, I n + 2, . . . , 2I n . The block tensor A B n is simply the concatenation of A and B at the 2n-th mode, which can be generalized to an arbitrary number of even-order paired tensors. More generally, given K even-order paired tensors X k ∈ R J ⊗I and a factorization K = K 1 K 2 . . . K N , the J ⊗ IK order mode row block tensor Y can be constructed in the following way: first, compute the 1-mode row block tensor concatenation over {X 1 , · · · , X K 1 }, {X K 1 +1 , · · · , X 2K 1 } and so on to obtain K 2 K 3 . . . K N block tensors denoted by X 1 , X 2 , . . . , X K 2 K 3 ...K N ; second, compute the 2-mode row block tensors concatenation over {X (1) 1 , · · · , X (1) 2K 2 } and so on to obtain K 3 K 4 . . . K N block tensors denoted by X 1 , X 2 , . . . , X K 3 K 4 ...K N ; third, keep repeating the process until the last N -mode row block tensor is obtained which we denote by Y = X 1 X 2 . . . X K , where IK = {I 1 K 1 , I 2 K 2 , . . . , I N K N }.
Analogously, one can also define the notions of n-mode column block tensors and mode column block tensors. Mode row/column block tensors possess many useful matrix-like properties such as the block tensor Einstein product which are useful in MLTI systems theory [9]. Moreover, the blocks of mode row/column block tensors map to contiguous blocks under ψ up to some permutations [9,51]. One can also define the inverse operation of extracting the component tensors from a given mode row block tensor. Let Y ∈ R J ⊗IK be a mode row block tensor constructed using the factorization K = K 1 K 2 . . . K N , then the component tensors X k ∈ R J ⊗I , k = 1, 2, . . . , K, can be extracted as: where, k = k 1 + N i=2 (k i − 1) i−1 l=1 K l for k n = 1, 2, . . . , K n . The colon operation : in (2.4) is a way to refer to the slices of a tensor as used in MATLAB.

Tensor decompositions.
There are a variety of notions of tensor decompositions such as CANDECOMP/PARAFAC decomposition, higher-order singular value decomposition, Tucker decomposition, tensor train decomposition and tensor singular value decomposition, which all play an important role in tensor algebra [7,9,12,16,27,28,43,46,47]. Of particular interest in this paper are the tensor train decomposition (TTD) and the tensor singular value decomposition (TSVD).
The TTD of an N -th order tensor X ∈ R J is given by where, • is the outer product, R = {R 0 , R 1 , . . . , R N } is the set of TT-ranks with R 0 = R N = 1, and X (n) ∈ R R n−1 ×Jn×Rn are called the core tensors of the TTD. Oseledets [43] showed that there exist optimal TT-ranks such that R n = rank reshape(X, reshape(X (n) , R n−1 J n , R n ) and X (n) = reshape(X (n) , R n−1 , J n R n ) are the left-and rightunfoldings of the core tensor, respectively [26]. Here I denotes the identity matrix, and reshape refers to the reshape operation in MATLAB. TTD has proved to be advantageous in both storage consumption and numerical stability, and truncating the TT-ranks results in a quasi-optimal approximation of X [43]. One can also define tensor trains for even-order paired tensors. Given an even-order paired tensor A ∈ R J ⊗I , the generalized TTD (GTTD) of A, also known as the TT-operator and TT-matrix defined in [16,44] respectively, is defined by where, A (n) ∈ R R n−1 ×Jn×In×Rn , and R is the set of GTT-ranks with R 0 = R N = 1 [9]. Clearly, TTD is a special case of GTTD with I = 1. One may even quantize the tensor trains (2.6) at each dimension, i.e., J n = J n1 J n2 . . . J nmn for some positive integer m n , in which a typical choice of J nmn is 2, in order to further reduce the complexity of TTD/GTTD and accelerate computations if the QTT-ranks are small [25,40,41]. We will refer to it as quantized TTD (QTTD). In fact, the QTTD of A ∈ R J ⊗I is the GTTD of the reshaped tensor A ∈ R (J 1 ⊗I 1 )×···×(J N ⊗I N ) where J n = {J n1 , J n2 , . . . , J nmn } and I n = {I n1 , I n2 , . . . , I nmn } such that |J n | = J n and |I n | = I n , respectively. A detailed algorithm for the conversion can be found in Algorithm A.1 from Appendix A. TSVD was first proposed by Brazell et al. [7] based on the isomorphism property and was extended by Sun et al. [56] for general non-paired even-order tensors. The results can be easily extended to even-order paired tensors. The economy-size TSVD (ETSVD) of an even-order paired tensor A ∈ R J ⊗I can be written as where, the superscript denotes the U-transpose operation, U ∈ R J ⊗R and V ∈ R I⊗R such that U * U = I and V * V = I (I denotes the U-identity tensor), S ∈ R R⊗R is an U-diagonal tensor containing the singular values of A along its diagonal S r⊗r , and R = R 1 R 2 . . . R N is equal to the unfolding rank of A. The ETSVD (2.7) also can be rewritten as where, X r ∈ R J ⊗1 , Y r ∈ R I⊗1 are the component tensors from the mode row block tensors U and V based on (2.4) with K n = R n for both tensors, and can be viewed as the left-and rightsingular tensors corresponding to the singular value σ r . All the singular values are arranged in descending order.

TT-algebra.
Computing the tensor algebraic notions related to even-order paired tensors using the unfolding ψ and matrix operations can be computationally demanding, especially when the size of tensors is large. Oseledets [43] showed that the basic linear algebra, such as addition, matrix-by-vector product and norms, can be done in the TT/GTT/QTTformat, without needing to go to the full tensor representation. The results were also extended for solving system of linear equations [18,42,44], and computation of Moore-Penrose (MP) inverse of unfolding matrices [26]. For simplicity, we summarize the key even-order paired tensor computations in the TT/GTT-format, and all the results also hold for the QTT-format.
Suppose that an even-order paired tensor A ∈ R J ⊗I is given in the GTT-format with cores A (n) and GTT-ranks R, the U-transpose of A can be obtained by transposing each component, i.e., We refer to an even-order "square" tensor A as weakly symmetric if A = A . Given two evenorder paired tensors A, B ∈ R J ⊗I in the GTT-format with cores A (n) , B (n) and GTT-ranks R, S, respectively, the element-wise TT-summation is given by A r n−1 ::rn , t n−1 = 1, 2, . . . , R n−1 , t n = 1, 2, . . . , R n B s n−1 ::sn , t n−1 = R n−1 + 1, R n−1 + 2, . . . , R n−1 + S n−1 , t n = R n + 1, R n + 2, . . . , R n + S n O, otherwise with T n = R n + S n for n = 1, 2, . . . , N − 1 and T 0 = T N = 1. Here O denotes the zero matrices with appropriate sizes. Given two even-order paired tensors A ∈ R J ⊗K and B ∈ R K⊗I in the GTT-format with cores A (n) , B (n) and GTT-ranks R, S, respectively, the TT-Einstein product is given by s n−1 ::sn ∈ R Jn×In , and t n = r n + (s n − 1)R n with T n = R n S n . The computational and memory complexities of the TT-Einstein product (2.11) are estimated as O(N J 3 R 4 ) and O(N J 2 R 4 ), respectively assuming J n = I n = K n ∼ J and R = S ∼ R where R can be viewed as the effective rank of GTTD defined in [45]. Furthermore, the Density Matrix Renormalization Group (DMRG) based algorithms proposed in [42,44] even enable one to solve the multilinear systems where, A ∈ R J ⊗J in the GTT-format and B ∈ R J in the TT-format. Readers may refer to [42,44] for more details.
If the GTT-ranks of even-order paired tensors are low, all the above computations can be achieved efficiently with low memory costs. Although the TT-summation and the TT-Einstein product between two even-order paired tensors may result in a new train with GTT-ranks larger than the optimal ones, one can apply TT-rounding (see Algorithm 2 in [43]) to resolve the problem with computational complexity O(N JR 3 ). In the following, we present some variations for computing mode row/column block tensors, and ETSVD for even-order paired tensors in the GTT-format. Again, the results also hold for the QTT-format if one first quantizes the tensor trains.

Block TT-format.
The block TT-format facilitates numerical computations for large scale problems [14,33]. Given two even-order paired tensors A and B of the same size in the GTT-format, we can construct the n-mode row block tensor without converting A and B to the full representation by filling zeros into the n-th cores. The main steps are summarized in Algorithm 2.1.
Algorithm 2.1 n-mode row block tensors 1: Given A, B ∈ R J ⊗I in the GTT-format with cores A (n) , B (n) and GTT-ranks R, S respectively, an integer n 2: Let A (n) = reshape(A (n) , R n−1 J n I n , R n ), B (n) = reshape(B (n) , S n−1 J n I n , S n ) The algorithm can be generalized to multiple blocks but would be more expensive in both computation and memory for a large number of blocks compared to the full format block tensor construction. This is because TT-summation would keep increasing the GTT-ranks requiring more memory space to store the cores. Even though one can apply TT-rounding during each TT-summation, the total computational cost of the algorithm is still very high for a large number of blocks. In addition, different ways of blocking may return different computational times and memory storages depending on the structure of GTTD. For example, the TTtoolbox function horzcat introduces extra modes to build block tensors for the TT-matrix class [45]. Therefore, careful choice of blocking algorithms can accelerate computations and save memory. On the other hand, the inverse operation of n-mode row block tensor also can be achieved in the GTT-format, i.e., the blocks of an n-mode row block tensor can be obtained by blocking the unfolding matrix of the n-th core, see Algorithm A.2 from Appendix A. Both Algorithm 2.1 and Algorithm A.2 can be extended to mode row/column block tensors given any factorization, and one can even choose the best factorization such that the mode row/column block tensor has the smallest GTT-ranks to achieve lowest memory requirements.
2.3.2. TSVD. Before discussing TSVD, we first introduce the notion of non-paired TTformat (NPTT-format). Given an even-order paired tensor A ∈ R J ⊗I in the GTT-format, its NPTT-format, naturally containing the information of its unfolding rank [9], is defined to be the TTD of the permuted tensor A np ∈ R J ×I . In the full representation, where, permute is the MATLAB permutation function. Significantly, the NPTTD of A can be constructed by manipulating the cores A (n) without converting back to the full format [45]. Detailed algorithms of this conversion are given in Appendix A, in which Algorithm A.3 transforms the GTTD of an even-order paired tensor to the TTD, and Algorithm A.4 permutes the TTD to the NPTTD. Note that both algorithms allow singular value truncations during the matrix SVD, which can remove redundancies and reduce the TT-ranks.
Klus et al. [26] exploited TTD to efficiently calculate the MP inverse of the matrix obtained from any chosen unfolding of a given tensor. By exploiting the relationship between the ETSVD and the unfolding matrix MP inverse, we adapt the framework of Klus et al. for the computation of the ETSVD of an even-order paired tensor, see Algorithm 2.2. In step 3, we can introduce a truncation threshold to obtain a low unfolding rank approximation of A, i.e., A K = K r=1 σ r X r * Y r . Appropriate truncations in the NPTT-conversion will return a good approximation of A. If such errors are negligible, A K can be regarded as optimal. After obtaining the left-and right-singular tensors, we can choose a factorization of R, such that the total memory cost of the GTTD is the smallest possible in order to construct the mode row block tensors U and V in (2.7). For computational convenience, we also prefer to set R N = R (or R 1 = R) in the factorizations.
The algorithm can be used to compute positive (or non-negative) U-eigenvalues and Ueigentensors (see details in [9]) for U-positive definite (or semidefinite) weakly symmetric tensors. In this case, all the singular values are U-eigenvalues, and the left-and right-singular tensors are equal to each other, i.e., X r = Y r , and are the U-eigentensors corresponding to the U-eigenvalues σ r . We will use this property to compute the U-eigenvalues and U-eigentensors of gramians, see subsection 3.1.2. Alternatively, the one-site optimization of the block Rayleigh quotient in the block TT-format proposed in [14] can also be used.
Instead of reshaping and computing the matrix SVD of ψ(A) ∈ R |J |×|I| which has at least O(min {|J | 2 |I|, |I| 2 |J |}) computational complexity, Algorithm 2.2 only computes a series of QR decompositions and matrix SVD's of the smaller-size left-and right-unfolding matrices of the cores. Hence, the computational complexity of the algorithm highly depends on the GTT-ranks and structure of A. Assume that J n = I n ∼ J and the GTT-rank R g ∼ R. If R remains unchanged or decreases during the NPTT-conversion with appropriate truncations, the computational and memory complexities are estimated by at most O(N 2 J 3 R 3 ) and O(N JR 2 ), respectively. Otherwise, both complexities may increase exponentially with N . As we will see in sections 3 and 4, dealing with large sparse MLTI/LTI systems, computing ETSVD through TTD could offer significant computational and memory benefits.
Algorithm 2.2 Economy-size TSVD 1: Given A ∈ R J ⊗I in the GTT-format with cores A (n) and GTT-ranks R g 2: Convert the GTTD of A to its NPTTD with cores X (n) and TT-ranks R using Algorithms A.3 and A.4 3: Compute the economy-size matrix SVD ofX (N ) , i.e.,X (N ) = UΣV for R = rank(X (N ) ), The left-and right-singular tensors of A are given by where, r − = {r 0 , r 1 , . . . , r N −1 } and r + = {r N +1 , r N +2 , . . . , r 2N } 6: return Left-and right-singular tensors U r and V r in the TT-format with singular values 3. MLIT system reduction/identification. The following MLTI system representation was first proposed by Chen et al. [8], where, X t ∈ R J is the latent state space tensor, Y t ∈ R I is the output tensor and U t ∈ R K is an input/control tensor. A ∈ R J ⊗J , B ∈ R J ⊗K and C ∈ R I⊗J are even-order paired tensors. Clearly, the MLTI system (3.1) can be transformed into an equivalent LTI system via ψ. The transfer function G(z) of (3.1) is defined by G(z) = C * (zI − A) −1 * B where z is a complex variable, and the superscript −1 denotes the U-inverse operation [9]. Based on (3.1), MLTI systems theoretic concepts including internal stability, reachability and observability are formulated in tensor forms. We recall several important notions from [8].
The MLTI system (3.1) is said to be reachable on [t 0 , t 1 ] if, given any initial condition X 0 and any final state X 1 , there exists a sequence of inputs U t that steers the state of the system from X t 0 = X 0 to X t 1 = X 1 . The pair (A, B) is reachable on [t 0 , t 1 ] if and only if the reachability gramian which is a weakly symmetric even-order square tensor, is U-positive definite. Here A t = A * A * t · · · * A. The infinite horizon reachability gramian can be computed from the tensor Lyapunov equation defined by The results of observability can be simply obtained by the duality principle. Chen et al. [9] also propose another equivalent MLTI representation with fewer parameters by using GTTD for efficient representation and computations, which is given by (we omit two colons in the cores) Moreover, one may apply QTTD to further reduce the number of parameters and complexity of MLTI systems if the QTT-ranks are small. We will see that QTTD is significantly advantageous in both computation and memory when dealing with large sparse systems/datasets. In many applications, such as high-dimensional partial differential equation problems, the states are high-dimensional tensors, while the number of inputs and outputs usually is much smaller than the number of states. It is thus desirable to approximate the large scale MLTI system with a lower-dimensional MLTI system. In this section, we propose three extensions of MLTI systems model reduction and identification methods -Higher-Order Balanced Truncation (HOBT), Higher-Order Balanced Proper Orthogonal Decomposition (HOBPOD) and Higher-Order Eigensystem Realization Algorithm (HOERA). Higher-Order Empirical Gramians (HOEG) will also be discussed. As mentioned, these model reduction and identification extensions also work for LTI systems if one first tensorizes the systems appropriately.
3.1. Higher-order balanced truncation. In LTI systems, balanced truncation (BT), introduced by Moore [37], is one of the most effective model reduction methods for stable linear input-output systems. BT is able to preserve system properties such as stability and passivity [48], but the exact balancing is expensive to implement for large systems requiring O(n 3 ) computational complexity and O(n 2 ) storage consumption [17]. We generalize BT to MLTI systems with fast computation and low storage consumption by exploiting the structure of the MLTI system (3.4). Like BT, the essence of HOBT is to find a higher-order transformation P ∈ R J ⊗J such that where, S ∈ R J ⊗J is the (unique up to permutation) U-diagonal tensor containing the Hankel singular values of the system which are independent of transformations based on the unfolding properties. We are first required to solve the tensor Lyapunov equations to obtain the reachability and observability gramians.
3.1.1. Solving Lyapunov equations. Solving Lyapunov equations by transforming them into large scale linear systems can be demanding both computationally and in terms of memory use. Many iterative methods for solving large scale Lyapunov equations, such as the Smith method and alternating direction implicit (ADI) iteration methods, are discussed in [55,60]. Nip et al. [39] proposed a novel approach by representing the large scale linear systems in the QTT-format, and solving the systems by the highly efficient DMRG based algorithms. We adapt here the Nip et al. approach to solving the algebraic tensor Lyapunov equation where, A, B ∈ R J ⊗J in the GTT/QTT-format, via the DMRG based solvers.
The key idea is to convert (3.6) into the following multilinear system where, the outer product between two GTTD can be readily obtained by simply connecting the two trains at the last and first cores, and B np is converted from B in the TT-format. The GTTranks of A • A are two copies of the GTT-ranks of A, and if B is sparse, the NPTT-conversion is fast. One of the advantages of this approach is that we can obtain the NPTT-format of the solution, i.e., X np , which can skip the NPTT-conversion in ETSVD for finding the Hankel tensor. Of course, the exact solution X ∈ R J ⊗J can be recovered from X np as desired in the GTT-format. Algorithm 3.1 applies the DMRG based linear TT-solver amen solver2 [45] to solve multilinear systems, which has proven highly efficient for low GTT/QTT-ranks inputs. However, there is no rigorous theoretical convergence analysis for DMRG based algorithms [44], so we will rely on numerical simulations to compare its performance.
where I is the U-identity tensor 4: Apply amen solver2 to solve the multilinear system L * X np = B np 5: return X np in the TT-format.

Computing Hankel singular values.
In LTI systems, after obtaining the reachability and observability gramians, computing the Hankel matrix requires computing the Schur decompositions of the gramians, which results in O(n 3 ) arithmetic operations [5]. It is therefore not suitable for large scale problems. Many low-rank methods such as Krylov subspace and ADI iteration have been proposed in [20,34] to resolve the problem. In the following, we present a fast and accurate algorithm to compute the Hankel singular values using ETSVD.
The reachability and observability gramians are U-positive semidefinite weakly symmetric tensors, so computing ETSVD of the two will return the U-eigenvalues and U-eigentensors of the gramians. Suppose that we obtain (W r ) np and (W o ) np from Algorithm 3.1. After left-and right-orthonormalizations, we can obtain the Cholesky-like factors of the two gramians, i.e., where, Z r ∈ R J ⊗C and Z o ∈ R J ⊗O are the mode row block tensors in the GTT-format for |C| ≤ |J | and |O| ≤ |J | consisting of the U-eigentensors multiplied by the square root of the corresponding U-eigenvalues, respectively. The computational complexity is estimated as O(N JR 3 ) assuming J n ∼ J and the TT-ranks of (W r ) np , (W o ) np ∼ R. In fact, if the original MLTI system possesses low GTT/QTT-ranks, |C|, |O| |J |. We therefore can define the Hankel tensor in the full format, i.e., where, ζ is the function that recovers a TTD/GTTD/QTTD in the full format. Then the Hankel singular values can be found more efficiently by reshaping and matrix SVD instead of using Algorithm 2.2. After calculating the Hankel singular values, HOBT can be achieved as described in Algorithm 3.2. In step 4, the left-and right-orthonormalization algorithms can be found in [26]. In step 5 and 7, one can choose the factorizations C n = |C|, O n = |O|, S n = |S| in C, O, S, respectively to return a LTI representation A r ∈ R |S|×|S| , B r ∈ R |S×|K| and C r ∈ R |I|×|S| . Since the system is in reduced form, it may not be necessary to represent in the GTT/QTT-format.
The reduced model is given by A r = Q * A * P, B r = Q * B and C r = C * P 10: return Reduced MLTI system A r ∈ R S⊗S , B r ∈ R S⊗K and C r ∈ R I⊗S .

Error bounds.
The regular BT provides an a priori error bound for the reduced system based on the Hankel singular values [15]. However, this error bound does not hold exactly for HOBT because several truncation errors occur during the GTTD/QTTD, the DMRG solver and the ETSVD. We nevertheless can estimate an error bound for HOBT based on the exact balancing. Suppose G and G r are the transfer functions of the full MLTI and reduced MLTI systems, respectively. Based on the triangular inequality, we have where, σ r are the exact Hankel singular values, S is the ETSVD truncation point, and is the H-∞ norm error between the two reduced systems from the exact balancing and HOBT.
Knowing that computing the exact error between G and G r can be computationally intense for large scale systems, (3.10) can be used to obtain an estimate of the performance of HOBT more efficiently. HOBT depends on the solutions of the tensor Lyapunov equations, but the DMRG based solvers are very sensitive to the GTT/QTT-ranks and the structure of MLTI systems. Hence, more efficient methods for approximating gramians are required.

3.2.
Higher-order empirical gramians. Instead of computing the gramians by solving the tensor Lyapunov equation (3.6) using Algorithm 3.1, one may directly compute them from data or numerical simulations [54]. This was the original approach used in Lall et al. [31,32] to extend balanced truncation to nonlinear systems. To obtain the higher-order empirical reachability gramain of a MLTI system with input size K, one runs the "higherorder impulse response" simulations of the primal MLTI system (3.1). Specifically, writing B = B 1 B 2 . . . B |K| based on the inverse operation of mode row block tensor construction (2.4), one constructs the state responses over t = 0, 1, . . . , T as follows: Consider T +1 = T 1 T 2 . . . T N snapshots of states for i = 1, 2, . . . , |K|, and arrange the snapshots in the following form: Then the higher-order empirical reachability gramian is given by In fact, this gramian is consistent with the definition in LTI systems. In other words, the unfolding ψ(W r ) is the empirical reachability gramian of the LTI system ψ(A), ψ(B), ψ(C) . According to Proposition 2.2 in [8], it follows that ψ( for a column permutation matrix P. Hence, ψ(W c ) ≈ ψ(X)ψ(X) = XPP X = XX . On the other hand, the procedure proceeds similarly for constructing the output snapshot tensor by collecting output snapshots Y ∈ R J ⊗IL from the simulations of the adjoint MLTI system (3.14)X t+1 = A * X t + C * V t , over t = 0, 1, . . . , L such that L + 1 = L 1 L 2 . . . L N .
Suppose we are only given the MLTI system (3.1) to begin with. After converting (3.1) to (3.4) by GTTD/QTTD with low GTT/QTT-ranks, all the steps described above can be efficiently achieved using TT-Einstein product and block TT-format for small values of T, L. However, if the impulse state responses data is provided for large T and L, it may be more efficient to construct X and Y using mode row block tensor in the full format.
3.2.1. Higher-order balanced POD. The method of balanced proper orthogonal decomposition (BPOD) was first proposed by Rowley [54] to deal with model reduction problems in high-dimensional input/output spaces such as in the context of fluid mechanics. Analogously, the goal of HOBPOD is to obtain an approximate HOBT that is computationally tractable for a large scale MLTI system. The method includes computing the balancing transformation from the collections of snapshots and a higher-order output projection method using ETSVD.
Suppose that the reachability and observability snapshots X and Y are given as in HOEG. Then the generalized Hankel tensor is defined by We can then simply apply Algorithm 2.2 to obtain the generalized Hankel singular values, and the remaining steps proceed similarly as HOBT. In Algorithm 3.3, we provide the case where the higher-order impulse state response data is given. In step 2, we use mode row block tensor to construct X and Y, but other choices of tensor blocking can also be used (e.g., extending the modes to build block tensors). Different choices of blocking would result in different compute times and memory storage requirements in the GTTD and ETSVD, but if one prefers to use QTTD, the choice of tensor blocking approach may not result in a significant difference. In particular, one can even use block matrices to form the generalized Hankel matrix H using ψ and apply QTTD directly to H. In step 4, one can also permute H to H np first in the full format and then apply QTTD to H np , but this approach may lose the flexibility of singular value truncations during the NPTT-conversion. In fact, it may be slower in computing the generalized Hankel singular values for similar level of accuracy. Additionally, depending on the situation, it may not be efficient to use the TT-Einstein product in step 4, 6, 7 to compute the generalized Hankel tensor H, transformations P,Q and the reduced system A r , B r , C r . For instance, in very high-dimensional MLTI systems, the GTT/QTT-ranks could be relatively large even if they are much less than the full ones, and so the TT-Einstein product may become a computational bottleneck. HOBPOD requires access to the MLIT system representation and its adjoint, and thus is not applicable in experimental setting. Moreover, the construction of the generalized Hankel tensor from the Einstein product between the reachability and observability snapshots X and Y can be computationally costly for high-dimensional MLTI systems. The eigensystem realization algorithm based approach alleviates these limitations of HOBPOD/BPOD.

Algorithm 3.3 Higher-order balanced POD
1: Given the snapshots of state responses X t ∈ R J andX l ∈ R J , two integer T, L and the MLTI system (3.1) A, B, C in the full format 2: Construct the reachability snapshot tensor using mode row block tensor where, X (k) = X 0 X 1 . . . X T and T + 1 = |T | 3: Construct the observability snapshot tensor Y ∈ R J ⊗IL using the adjoint state responseŝ X l similarly as step 2 4: Compute the generalized Hankel tensor H = Y * X ∈ R IL⊗KT and apply GTTD (or QTTD) to H 5: Compute the unfolding rank S approximation of H using Algorithm 2.2, i.e., The reduced model is given by A r = Q * A * P, B r = Q * B and C r = C * P 8: return Reduced MLTI system A r ∈ R S⊗S , B r ∈ R S⊗K and C r ∈ R I⊗S .

3.3.
Higher-order ERA. The eigensystem realization algorithm (ERA) was first proposed in [24] as a model identification and reduction tool for LTI systems. It was later shown [36] that ERA is theoretically equivalent to BPOD for discrete time systems but with significantly lower computational cost. Similarly, the idea of HOERA is to construct a generalized Hankel tensor using higher-order impulse response simulations or experiments without having access to the MLTI system (3.1).
First, we need to collect the snapshots of the outputs Y k ∈ R I⊗K and form the generalized Hankel tensor as follows: where, Y k = C * A k * B are called Markov parameters, T +L+1 is the number of snapshots with T + 1 = T 1 T 2 . . . T N and L + 1 = L 1 L 2 . . . L N . T and L are usually chosen to be sufficiently large. By applying Algorithm 2.2, we can obtain the generalized Hankel singular values, and the remaining steps are summarized in Algorithm 3.4. Like HOBPOD, the QTTD of H is invariant under different choices of tensor blocking. We can also convert the tensors to the full format after ETSVD in computing the reduced model at step 4 in case the TT-Einstein product becomes a computational bottleneck. The reduced MLTI system obtained from HOERA is consistent with that obtained from ERA on the unfolding LTI system. In other words, the reduced system ψ(A r ), ψ(B r ), ψ(C r ) can be achieved from ERA on the LTI system ψ(A), ψ(B), ψ(C) . Due to Proposition 2.2 in [8], it follows that where, Q is a row permutation matrix and P is a column permutation matrix. Suppose that the matrix SVD on H is given by H = USV . Then ψ(U) = QU, ψ(V) = P V and ψ(S) = S. Hence, where, ψ(H 1 ) = QH 1 P, and Col first (H) and Row first (H) represent the first block column and row of the generalized Hankel matrix H. Moreover, if the outputs Y k are collected from "higher-order step response" simulations or experiments, HOERA can still be applied. It can be shown that the difference between two consecutive outputs is given by where, B i is the i-th row tensor block of B activated by the control inputs for 1 ≤ i ≤ |K|. Hence, one can apply Algorithm 3.4 to the snapshots {δY 1 , δY 2 , . . . , δY T +L+1 } and obtain the reduced MLTI model. Algorithm 3.4 Higher-order eigensystem realization algorithm 1: Given the snapshots of the output Y k ∈ R I⊗K , two integers T, L 2: Construct the generalized Hankel tensor H ∈ R IL⊗KT using mode block tensor with T + 1 = |T | and L + 1 = |L| and apply GTTD (or QTTD) to H 3: Compute the unfolding rank S approximation of H using Algorithm 2.2, i.e., H ≈ U * S * V with U ∈ R IL⊗S , V ∈ R KT ⊗S , S ∈ R S⊗S where S = {S 1 , S 2 , . . . , S N } 4: The reduced model is given by in the full format where Col first (H) and Row first (H) represent the first block column and row of H and

Synthetic datasets.
In this example, we consider a multiple input multiple output (MIMO) MLTI system (3.4) with low QTT-ranks random tensors in the QTT-format A ∈ R 2⊗2 , B ∈ R 2⊗2 and C ∈ R 2⊗2 . Consequently, the number of states in the unfolding LTI corresponding to MLTI system (3.1) is 2 N . We compare the computing times of HOBT and ordinary BT for different values of N (and hence number of states) in Figure 1 (A), and the error bounds in Figure 1 (B). As can be seen, the computing time for ordinary BT grows rapidly with the number of states N but remains bounded for HOBT, while the errors in model reduction using the two approaches are similar. For the error bound, we use the sum of the residual Hankel singular values for the ordinary BT, and the bound (3.10) for HOBT.

2D heat equations with control.
The heat equation is a partial differential equation that describes the evolution of heat distribution in a solid medium over time. The 2D heat equation on the square D = [−π, π] 2 with localized point control and Dirichlet boundary conditions is given by where, c > 0, u t ∈ R is a one-dimensional control input and δ(x) is the Dirac delta function centered at zeros. In this example, we apply HOBPOD and HOERA to the discretized heat equation (4.1) to find the reduced balanced system.
We use a second-order central difference to approximate the Laplacian, first-order difference in time, and we approximate the Dirac delta as a Kronecker delta function at the nearest grid point with unit mass, leading to a MTLI system of the form (3.1), where X t ∈ R N ×N is the 2D-temperature field at instance t, and A ∈ R N ×N ×N ×N , B ∈ R N ×1×N ×1 and C ∈ R 1×N ×1×N are system tensors. The tensor A is the tensorization of matrix c 2 ∆t h 2 ∆ dd ∈ R N 2 ×N 2 , where ∆ dd is the discrete Laplacian on a rectangular grid with Dirichlet boundary conditions, and B is the tensorization of 1 h 2δ (x), whereδ(x) is a vector which is zero everywhere except at the entry corresponding to the grid point closest to the origin (see [39] for details). Here, h = ∆x = ∆y are spatial resolution of grid sizes in the x and y directions, respectively, ∆t is discretization in time, and we ensure the Courant-Friedrichs-Lewy (CFL) condition c 2 ∆t h 2 < 1 is met for numerical stability.
We assume measurement at a single discrete location, i.e., I 1 = I 2 = 1, and generate forward/adjoint input/output snapshot data over the T + L + 1 time steps for different values of spatial discretization δx. Let N = 2π/δx be the dimensions of the two modes of the temperature tensors X t . Figure 2 (A) shows the comparison of HOBPOD/HOERA computational times with BPOD/ERA as a function of system dimension, i.e., number of states. We consider both GTTD and QTTD approaches for HOBPOD/HOERA, and both are found to be more efficient than BPOD/ERA. Furthermore, as expected, QTTD is more efficient than GTTD if the tensor train conversion time is included. Also shown in the figure are the curves for GTTDw (and similarly QTTDw), which stands for computational time for GTTD based HOBPOD/HOERA without accounting for time to convert the tensors to their corresponding GTT-format. This is because some special tensors/matrices like the Laplacian can be efficiently constructed in the GTT/QTT-format (e.g., see the function tt qlaplace dd in the MATLAB TT-Toolbox [45]). Note that while this conversion time can be significant for GTTD, for QTTD it is found to be negligible. Figure 2 (B) shows difference between transfer functions obtained via BPOD/ERA and HOBPOD/HOERA, indicating that HOBPOD/HO-ERA obtains the reduced model with similar accuracy to BPOD/ERA.

Cancer cell image video dataset.
The Fas-Associated Death Domain (FADD) is an adaptor protein known for its role in cell extrinsic apoptosis (cell death), and studies have revealed that deletion of FADD is able to suppress lung cancer development [6]. In this example, we try to capture the cell cycle dynamics of a cancer cell with and without downregulation/suppression of FADD by HOERA/ERA. The downregulation/suppression of FADD is accomplished by the intranasal inhalation of an adenovirus expressing Cre recombinase (AdCre) [6], which can be viewed as a step input to the cell cycle dynamics.
The dataset contains two cancer cell image videos, see Figure 3. One is a normal cancer cell, and one is a cancer cell with deletion of FADD. Each video has 145 frames with pixels 64 × 64 and two channels. The channels are used to indicate cell cycle state transitions. For both videos, we binarize each frame and construct the generalized Hankel matrix H ∈ R 2 17 ×128 with T = 127 and L = 15 using ψ. Then we apply QTTD to H and obtain a tensor train with J n = 2 for n = 1, 2, . . . , 17, I n = 1 for n = 1, 2, . . . , 10 and I n = 2 for n = 11, 12, . . . , 17. In the following study, we compare the Hankel singular values in HOERA and ERA with different truncations in the QTTD step and observe the dominant features from the first few left-singular tensors and left-singular vectors, respectively, for the two videos. Along the way, we report the number of parameters that require to store the outputs of ETSVD and economy-size matrix SVD. The results are shown in Figure 4. The ETSVD has a huge advantage in memory saving compared to the economy-size matrix SVD, even if we recover some left-singular tensors in the full format in order to show the dominant features. The Hankel singular values have a similar decay pattern to those computed from matrix SVD for smaller truncation thresholds . The images of the first few left-singular tensors from the two videos quantitatively show that the protein FADD can keep the cancer cell stationary during the cell cycle.

4.4.
Room impulse responses dataset. The reduction of room reverberation has become increasingly significant in human life. In this example, we try to apply HOERA/ERA to a real world database of binaural room impulse responses (BRIR), referred to as the Aachen Impulse Response (AIR) database [22,23]. The goal is using HOERA to efficiently and accurately capture the acoustic dynamics of a given room, which may potentially help promote the development of algorithms for dereverberation. The BRIR is a one-dimensional output signal, i.e., scalar, so one can simply apply ordinary ERA to the dataset. However, we observe that the generalized Hankel matrix built from the impulse response is huge given the large sampling frequency, and it also contains many extremely small values from the later time points. These suggest that HOERA may be more suitable in estimating the acoustic dynamics.
The AIR database includes many diverse scenarios. In the following studies, we focus the attention on the case where the room impulse responses are collected by a bottom microphone with a bottom-top mock-up phone in an office (see details in [22,23]). We compare the computational times for identifying the acoustic systems using HOERA and ERA with certain numeric accuracies. Let the number of snapshots be T = L = 2 N − 1 for a positive integer N . Then we apply QTTD to the generalized Hankel matrix H ∈ R 2 N ×2 N and obtain a tensor train with J n = I n = 2 for n = 1, 2, . . . , N . The remaining steps are followed as discussed in Algorithm 3.4 to find the approximating acoustic model.
The results are given in Tables 1 and 2. In the first study, with the increase of T and L, we can allow larger prescribed truncation thresholds in the NPTT-conversion (Algorithm A.4), which accelerates the computation of ETSVD. HOERA exhibits computing time advantage in estimating the system model with low relative errors for N ≥ 12, see Figure 5. Of course, ERA can provide better relative errors with the same number of singular values retained as in HOERA. In the second study, when N = 13, larger prescribed truncation thresholds in the NPTT-conversion enable faster computations and lead to similar relative errors compared to ERA, see Table 2. Note that we report the times of computing ETSVD with the times of QTTD and economy-size matrix SVD for HOERA and ERA, respectively, in the both studies. Table 1 Computational times of HOERA and ERA for different values of N with sampling frequency 5600. "Threshold" represents the maximum truncation threshold required in the NPTT-conversion to maintain the relative errors less than 0.1. "Relative Error 1 " and "Relative Error 2 " represent the relative errors between the real data and the impulse responses generated by the models identified by HOERA and ERA, respectively. 5. Discussion. The numerical studies reported in the previous section highlight the computational efficiency of different flavors of the proposed HOBT approaches compared to the analogous standard matrix based BT for model reduction/identification of high-dimensional MLTI/LTI systems. These gains come from low rank/sparse structure in the data/underlying system which tensor decomposition methods such as the TTD exploit. While TTrepresentation and associated algebra are at the core of our proposed framework, it may not always be beneficial to work with the TT-format. In particular, as pointed out in subsections 3.2.1 and 3.3, if the number of snapshots, T or L, is large, computing the block TT-format and TT-Einstein product can become computational bottlenecks, and it may be preferable to work with the full tensor format or even the unfolded matrix representation. In some circumstances, it may still be possible to control growth of TT-ranks by appropriate truncation, choice of blocking operations and choice of appropriate factorizations used in the block tensor construction as discussed in subsection 3.2.1. From our numerical studies, we found that QTTD is preferable to TTD as it is less sensitive such choices. Additionally, as proposed in subsections 3.1.2, 3.2.1 and 3.3, one can switch to matrix representation/computation using unfolding at appropriate stages in the algorithm. However, more theoretical and numerical investigation is required to assess under what conditions TT-format versus full tensor format versus matrix representation is computationally more efficient, and is an important avenue of future research. Finally, it should be noted that while the proposed framework is naturally suited for MLTI systems, it can also be applied to the standard vector based LTI system representation by tensorizing it into a suitable tensor form or applying QTTD directly as we demonstrated in the room impulse response example.
The proposed data-guided TTD based BT methods can be applied to many large biological and engineering systems/datasets. For example, cellular reprogramming, which involves introducing transcription factors as a control mechanism to transform one cell type to another, requires a comprehensive view of genome-wide structure, gene expression, the proteome and phenotype fitted in a tensorial representation [52,57]. Furthermore, tensor analysis has played an increasingly significant role in machine learning recently [10,61]. The cancer cell image video example using HOERA/ERA is similar to the eigenfaces problem [59], which indicates that the framework of computing ETSVD using TTD would also aid in machine learning algorithms such as principal component analysis, regression analysis and spectral graph theory for large sparse datasets. In particular, ETSVD can help accelerate the computation of von Neumann entropy for large sparse graphs like aviation and citation networks for example.
6. Conclusion. In this paper, we developed a Higher-Order Balanced Truncation, Higher-Order Balanced Proper Orthogonal Decomposition and Higher-Order Eigensystem Realization Algorithm based model reduction and identification approach for input/output multilinear/linear time invariant systems. The proposed framework exploits tensor decompositions such as tensor train decomposition (including standard TTD, generalized TTD and quantized TTD) and economy-size tensor singular value decomposition to compress the tensors which facilitate efficient computations while retaining the accuracy. We applied our data-guided model reduction/identification framework to real world biological and engineering systems and achieved outstanding performances in computational efficiency and memory consumption with rigorously quantified errors. We are also currently developing a Dynamic Mode Decomposition (DMD) type variant of our framework. Some initial work generalizing DMD using TTD appears in [26]. It would also be worthwhile to develop an observer and feedback control design framework and associated tensor based computational techniques for estimation/control of MLTI/LTI systems and apply these techniques in real world complex systems. Nonlinearity and stochasticity in tensor based dynamical system representation and analysis framework are also important for future research.

6:
Compute economy-size SVD of A, i.e., A = USV with S ∈ R R q nk ×R q SetÃ (nmn) =Ã 11: end for 12: return A ∈ R J ⊗I in the QTT-format with coresÃ (nk) and QTT-ranks R q .
Algorithm A.2 Inverse n-mode row block tensors 1: Given C ∈ R J ⊗I in the GTT-format with cores C (n) and GTT-ranks R, an integer n such that I n is even 2: Set C (n) = reshape(C (n) , R n−1 J n I n , R n ) 3: Let A and B be the generalized tensor trains copied from C except the n-th cores A (n) = reshape(C I n , R n ), respectively 4: return A, B ∈ R J ⊗L in the GTT-format with cores A (n) , B (n) and GTT-ranks R, respectively where L = I except L n = 1 2 I n .
Algorithm A.3 Conversion of GTTD to TTD 1: Given A ∈ R J ⊗I in the GTT-format with cores A (n) and GTT-ranks R g 2: for n = 1, 2, . . . , N do 3: Set A = reshape(A (n) , R g n−1 J n , I n R g n )

4:
Compute economy-size SVD of A, i.e., A = USV with S ∈ R Rn×Rn