{'fraction_non_alphanumeric': 0.042796379543209626, 'fraction_numerical': 0.012840929707500958, 'mean_word_length': 3.917525773195876, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Orthogonal Decomposition of Tensor Trains

In this paper we study the problem of decomposing a given tensor into a tensor train such that the tensors at the vertices are orthogonally decomposable. When the tensor train has length two, and the orthogonally decomposable tensors at the two vertices are symmetric, we recover the decomposition by considering random linear combinations of slices. Furthermore, if the tensors at the vertices are symmetric and low-rank but not orthogonally decomposable, we show that a whitening procedure can transform the problem into the orthogonal case. When the tensor network has length three or more and the tensors at the vertices are symmetric and orthogonally decomposable, we provide an algorithm for recovering them subject to some rank conditions. Finally, in the case of tensor trains of length two in which the tensors at the vertices are orthogonally decomposable but not necessarily symmetric, we show that the decomposition problem reduces to the novel problem of decomposing a matrix into an orthogonal matrix multiplied by diagonal matrices on either side. We provide and compare two solutions, one based on Sinkhorn's theorem and one on Procrustes' algorithm. We conclude with a multitude of open problems in linear and multilinear algebra that arose in our study.


Introduction
With the emergence of big data, it is more and more often the case that information is recorded in the form of a tensor (or multi-dimensional array). The importance of being able to decompose such a tensor is (at least) twofold. First, finding the decomposition provides hidden information about the data at hand, and second, having a concise decomposition of the tensor allows us to store it much more efficiently. One of the biggest obstacles in dealing with tensors, however, is that decomposing them is often computationally hard. For example, finding (the number of terms of) the CP-decomposition [1] of a general tensor is NP-hard [2], and the set of tensors of CP rank at most r is not closed for any r ≥ 2, making the low-rank approximation problem impossible in some instances. Thus, it is sometimes hard to use this decomposition despite its widespread applicability.
While there exist algorithms for finding the CP decomposition of a general tensor, such as Alternating Least Squares [3], they can often get stuck at local optima and cannot guarantee convergence to a global optimum. Finding the CP-decomposition of special subclasses of tensors, however, can be done efficiently. For example, low-rank tensors can be decomposed via Jennrich's algorithm and their subfamily of orthogonally decomposable tensors can be decomposed via the slice method [4] or via the tensor power method [5]. Furthermore, the set of orthogonally decomposable tensors of bounded rank is closed [6], and their eigenvectors and singular vector tuples can be found efficiently [7,8], making the family of orthogonally decomposable tensors as appealing as the set matrices. However, such tensors are very rare -the rank of a general n × n × · · · × n (d-times) tensor is O(n d−1 ), while that of an orthogonally decomposable one is n. The recently developed tensor subspace power method [9] allows for efficient decompositions of symmetric tensors of rank up to O(n When the tensor train has length two (cf. Figure 4), and the orthogonally decomposable tensors at the two vertices of the network are symmetric, we show how to recover the decomposition by considering random linear combinations of slices (cf. Section 3). Furthermore, if the tensors at the vertices are symmetric but not orthogonally decomposable, we show that a whitening procedure can transform them into orthogonally decomposable ones yielding a way of decomposing the tensor. This type of tensor network decomposition is equivalent to a hidden variable graphical model decomposition, and our algorithm can be applied to recover the hidden information (cf. Section 3.4). When the tensor network has length three or more and the tensors at the vertices are symmetric and orthogonally decomposable, we provide an algorithm for recovering them as long as their ranks obey certain inequalities (cf. Section 4). Finally, in the case of tensor trains of length two in which the tensors at the vertices are orthogonally decomposable and not necessarily symmetric, we show that the decomposition problem reduces to the problem of decomposing a matrix into an orthogonal matrix multiplied by diagonal matrices on either side (cf. Section 5). In Section 6 we show how to find such a matrix decomposition in the full-rank case using Sinkhorn's theorem [14] or using Procrustes' algorithm [15], and show how the latter solution can be generalized to any rank case. This leads to an application in electronic structure theory for finding an orthogonal factorization of the electron repulsion integral (ERI) tensor, in the case when such a factorization exists (cf. Section 5.1). We report our experiments in Section 7. We conclude in Section 8 with a multitude of open problems in linear and multilinear algebra that arose in our study. The sections of this paper are organized according to the flowchart in Figure 1.

Background
In this section we provide the necessary background and notation on tensors, tensor decompositions, tensor networks, and orthogonally decomposable tensors.

Tensors
We denote the set {1, ..., n} as [n], for n ∈ N. A d-tensor with real entries is an element in R n1×···×n d for some n 1 , . . . , n d ∈ N. Scalars, vectors, and matrices are 0-, 1-, and 2-tensors, respectively. We denote scalars by lowercase letters (e.g. λ), vectors by bolded lowercase letters (e.g. u), matrices by bolded uppercase letters (e.g. U), and d-tensors by script uppercase letters (e.g. T ), for d ≥ 3. If u 1 ∈ R n1 , ..., u d ∈ R n d , then one can form their tensor product, which is the d-tensor T = u 1 ⊗ ... ⊗ u d ∈ R n1×...×n d with entries Such tensors are called rank-1 tensors. Every tensor is a finite linear combination of rank-1 tensors, and the smallest number of such terms in the combination is called the rank of the tensor.

Tensor Networks
There is an operation generalizing matrix multiplication to tensors: If T ∈ R n1×...×ni−1×ni×ni+1...×n d and S ∈ R m1×...×mj−1×mj ×mj+1...×m such that n i = m j , then their contraction along the i th and j th modes is the tensor R ∈ R n1×.  Figures 2a and 2b, where the dimension of each vector component labels an edge. If it is clear from the context, these edge labels can be dropped. Note that regardless of the number of tensor product terms in the sum defining T and S, the same diagram of a vertex with three edges is used to represent them; what the diagram illustrates is the dimensionality of T and S, as they are formed by tensor products of three vectors. If n c = n f , then T and S can be contracted along the c-and f -modes to form the tensor Figure 3, and n e = n , then we can contract R with the tensor 2c along the e-and -modes to form the tensor The tensor networks in Figures 2c and 3b are called tensor trains of length 2 and 3, respectively, and generalize to longer trains. Tensor trains are generated by the sequential contraction of 3-tensors, like T , Figure 3: Carriage R and length 2 tensor train 2c contract to form a length 3 tensor train S, and R. These 3-tensors at the vertices of the network will be referred to as carriages. Tensor trains can produce tensors of high rank: ij ) ⊗3 (j = 1, . . . , L) at its vertices, each of rank at most n, is at most n L = n d−2 , where d is the order of the tensor.
Proof. Indeed, such a decomposition yields the following tensor which has rank at most n L and order d.
Matrix product states often appear in the condensed matter physics literature with one additional vertex and edge attached to each end of the diagram, as opposed to our diagrams as in Figure 2 [10,11]. This format can be converted to the format introduced here simply with a contraction of the single edges, thus removing the vertices at each end. Conversely, given a tensor network in the form of our diagram, a multiplication by the identity matrix at each end will yield the matrix product state form. Thus, as the two network formats are equivalent, we proceed with the use of the term "tensor train" to describe the networks we study in this paper. The reason for the modified form is to better generalize and extend the class of tensors known as orthogonally decomposable tensors, as studied in [4,5,6,7,8].

Orthogonally Decomposable Tensors and Tensor Trains
i . Orthogonally decomposable tensors are appealing because they can be decomposed efficiently via the slice method [4] or via the tensor power method [5]. For more on orthogonally decomposable tensors, please refer to [5,6,7,8]. A tensor train is said to be symmetric/orthogonal if it is generated by carriages that are symmetric/orthogonally decomposable. The goal of this paper is to decompose such tensors.
⊂ R nc are each an orthonormal set. r T is said to be the rank of T ,â i ,b i ,ĉ i the vectors, and λ i the coefficients of T . To decompose a tensor train, we therefore mean to find the ranks, vectors, and coefficients of each carriage which generate it. Our goal is to decompose tensors which are known to have an orthogonal tensor train form; we are not aware of a procedure for testing whether a given tensor has an orthogonal tensor train decomposition. The methods we present in Sections 3, 4, and 5 for decomposing tensor trains can be used reliably when there is no noise in the tensor. In Section 7, we perform numerical experiments to test our methods both in the absence and presence of noise. One of the main ambitions of tensor network research is to find efficient decompositions. Tensor network diagrams can have line segments, which correspond to tensor trains and is A B Figure 4: A tensor network diagram of T what we study in this paper, as well as loops, which correspond to tensor rings and have been studied in [16]. For more on tensor networks, please refer to [10,11].
Lastly, there are objects which we will encounter, such as scalars and orthonormal sets of vectors, which lie in a Zariski-closed subset of R n , for some n pertaining to each object. We will say that such an object is generic if it lies in some Zariski-open subset of this closed set. In other words, if an object is randomly drawn from this closed set, then with probability 1, it will be drawn from the open subset. We also abbreviate "up to permutation and sign" as UTPS.

Symmetric Orthogonal Decomposition of Tensor Trains of Length 2
In this section we investigate the problem of decomposing an order-4 tensor T according to a tensor train network (cf. Figure 4) such that the carriages A and B are symmetric. We begin in Section 3.1 by also assuming that A and B are symmetric orthogonally decomposable, and we show in Section 3.2 how to use whitening in the non-orthogonal case.

Symmetric Orthogonal Decomposition of Length 2
Problem 1. Let T ∈ R n×n×n×n be a 4-tensor admitting the following decomposition ⊂ R n are generic orthonormal sets, and λ i , µ j ∈ R \ {0}. In other words, assume that T is a tensor train generated by Given T as in (1), find the above decomposition, including all ranks, vectors, and coefficients.
We solve this problem by adapting Kolda's slice method [4]. Consider generic weighted sums over the last two indices and first two indices of T for some generic α i3i4 , β i1i2 , σ i , γ j ∈ R. Here Σ, Γ ∈ R n×n are diagonal with σ i or 0 and γ j or 0 along their diagonals, respectively, and U, V ∈ R n×n are orthogonal matrices with u i and v j in r A and r B of their columns, respectively. The vectors u i and v j can, therefore, be found UTPS via an eigendecomposition on S A and S B , since they correspond to the nonzero eigenvectors. Thus, we have found all ranks and vectors. To find the coefficients, construct the matrix R ∈ R r A ×r B with entries Note that uî, vĵ = 0 since u i and v j are generic. Then, R = λµ is a rank-1 matrix, where λ and µ are the vectors whose components are λ i and µ j . Performing Singular Value Decomposition (SVD) on R, if ζ and η are the first left and right singular vectors and τ is the largest non-zero singular value, then λ = τ ζ and µ = η is a solution. This non-uniqueness is due to non-zero scaling.
Pseudocode for this procedure can be found in Algorithm 1.

Whitening: Symmetric Non-Orthogonal Decomposition of Length 2
If T has a decomposition as in (1) where {u i } r L i=1 and {v j } r R j=1 are linearly independent but not orthonormal sets, then we can adapt Kolda's method of whitening [4]. The generic weighted sums over indices of T have two different decompositions: one in terms of the vectors u i and v j , and one in terms of its "skinny" eigendecomposition: are diagonal with non-zero diagonal entries, U ∈ R n×r A and V ∈ R n×r B have u i and v j as columns, and X A ∈ R n×r A and X B ∈ R n×r B are orthogonal matrices. If C A and C B are positive semi-definite (PSD), then Σ, D A , Γ, and D B have positive diagonal entries. Hence, let Then we see, for instance, ×r A has orthogonal columns, and similarly for V = W B V ∈ R r B ×r B . This means that the tensor train obtained by contracting T with the matrix W A along its left two dangling edges and with the matrix W B along its right two dangling edges: for appropriate λ i , µ j ∈ R, is symmetric and orthogonal and can be decomposed by the method described earlier. If u i are the vectors found from the eigendecomposition, then one can recover A u i and similarly for v j . Algorithm 1 provides pseudocode for decomposing symmetric tensor trains of length 2 in either the orthogonal or non-orthogonal case, and its correctness is summarized by the following result: ⊂ R n are generic vectors, then Algorithm 1 recovers r A and r B , and λ i , Algorithm 1 Symmetric Decomposition of Tensor Trains of Length 2 Input: T ∈ R n×n×n×n that decomposes according to a symmetric tensor train of length 2. An indicator if A and B are non-orthogonal, and therefore we need to apply whitening first. Output: ρ, τ ← generic real n × n matrices 4:

the coefficients and (orthonormal) vectors forming
] 31: λ ← largest singular value multiplied by corresponding left singular vector from SVD of R 32: µ ← corresponding right singular vector from SVD of R

Remarks
It is not guaranteed that C A and C B are PSD, and hence if they are not, we recompute generic weighted sums over indices of T until they are PSD or exit with failure. We analyze the performance of Algorithm 1 in Section 7. It is also straightforward to generalize to the case where each carriage in the train is a symmetric m-tensor, for some m > 3. Furthermore, if there are p contracted edges between A and B, then one simply divides by u i , v j p when forming the matrix R from Section 3.1.

Applications
The tensor decomposition discussed in Section 3.2 is equivalent to the decomposition of the joint distribution of the four leaf variables in the graph in Figure 5 below (note that the joint distribution of four variables each taking n values is precisely an n × n × n × n tensor). This follows directly from the duality between tensor networks and graphical models in [13]. Note that the three hidden variables take r, s, and n values respectively, and finding the decomposition discussed in Section 3.2 recovers the joint distribution of all the variables in the graphical below. r s n hidden n n n n Figure 5: A hidden variable graphical model equivalent (in the sense that they yield the same tensor decomposition) to a length-2 tensor train whose carriages are symmetric and have ranks r and s, respectively.

Symmetric Orthogonal Decomposition of Tensor Trains of Length
L ≥ 3 In this section we consider longer tensor trains such that the carriages are symmetric and orthogonally decomposable. Figure 3: A tensor network diagram of T Problem 2. Let L ≥ 3 and T ∈ R n L+2 be an (L + 2)-tensor admitting the following decomposition ij =1 ⊂ R n is a generic orthonormal set, and λ j ij ∈ R are generic. In other words, assume that T (cf. Figure 4) is a tensor train generated by Given T , find the decomposition (2), including all ranks, vectors, and coefficients.
The assumption that {x j ij } r X j ij =1 and λ j ij are generic will be important for our results on how to decompose T , as the results may not hold otherwise. Because of the genericity assumption, our results hold with probability 1. We present a solution when T satisfies the following condition: The DRC essentially states that there is not too much information lost about the ranks, vectors, and coefficients of the carriages upon contraction. The steps of the solution involve a method called kernel completion, which we describe in detail later, and are summarized below. We present pseudocode in Algorithm 2.
where v ∈ R n is a generic vector. Then sequentially apply kernel completion starting from X 2 and ending at X L−1 , obtaining a collection of orthonormal vectors where v ∈ R n is a generic vector. Then sequentially apply kernel completion starting from X L−1 and ending at X 2 , obtaining a collection of orthonormal vectors {x L−1,RL based on whichever of r Xj ,LR or r Xj ,RL is greater. If they are equal, choose either set.
4. Construct an L-tensor R ∈ R r X 1 ×...×r X L whoseî 1 ...î L -entry is given by Apply the rank-1 alternating least squares [3] algorithm on R to obtain all the coefficients. We now explain each of these steps in more detail.
Step 1: Decomposing the train from left to right We can always find the vectors at the "ends" of the train {x 1 i1 } and {x L i L } UTPS regardless of whether or not T satisfies the DRC. Let v ∈ R n be a generic vector and consider X 1 = T (·, ·, v, ..., v) i1 } can therefore be found UTPS via an eigendecomposition of X 1 , as they correspond to the non-zero eigenvalues. The same method can be used to find X L ∈ R n×n , which has {x L i L } in r X L of its columns, UTPS. Note that {x 1 i1 } and {x L i L } can also be found using Kolda's method as with trains of length 2, which forms an alternative X 1 by constructing generic weighted sums over all but two indices of T where γ i3,...,i L+2 ∈ R are generic. We found our method to be efficient and simple to implement for long trains, hence we use it.

Kernel completion
Complete the set of orthonormal vectors {x 1 i1 } r X 1 i1=1 found above to an orthonormal basis and let X 1 ∈ R n×n be the orthogonal matrix whose columns are these basis vectors. Denote by X 2 ∈ R n×n an orthogonal matrix whose first r X2 columns are {x 2 i2 }, let v ∈ R n be a generic vector, and consider T = T (v, ·, ·, v, ..., v) Multiplying on the left by X 1 and on the right by X 1 , we get We seek to find a diagonal matrix L = diag 1 , ..., r X 1 , 1, ..., 1 ∈ R n×n with which to multiply on the left of the above equation so that the non-zero entries of Λ get canceled.
Since QMQ is symmetric, the resulting matrix will have a symmetric top-left r X1 × r X1 block. If r X1 = 1, this symmetrizing procedure is unnecessary, since a top-left 1 × 1 corner block of any matrix is symmetric. This condition allows us to find a linear system in the unknowns 1 , . . . , r X 1 which has the following form. Let T ∈ R r X 1 ×r X 1 be the top-left r X1 × r X1 corner block of the matrix X 1 TX 1 . Form the matrix L L L ∈ R ( r X 1 2 )×rX 1 , where for each row we choose two indices i, j ∈ [r X1 ], i = j, let the i th entry of the row equal T i,j , the j th entry of the row equal −T j,i , and the remaining entries be 0. Then we have the following result: Proposition 1. If T is generated by orthogonally decomposable carriages whose vectors and coefficients are generic, then nullsp(L L L ) = span((λ −1 1 , ..., λ −1 r X 1 )), which gives us precisely 1 , . . . , r X 1 .
Proof. Let R ∈ R r X 1 ×r X 1 be the top-left r X1 × r X1 corner block of QMQ , which is symmetric: Then by the definition of L L L , ( 1 , ..., r X 1 ) ∈ nullsp(L L L ) satisfies Since we assumed that all orthonormal sets and coefficients are generic, we have for some constant C ∈ R, we prove the following lemma: Proof. We proceed by induction on n. When n = 1, the condition 1 1 = 1 = λ1 λ1 is vacuous and there must exist a C ∈ R such that 1 = C λ1 . Now suppose the statement is true for a sets of size n. Given satisfying the conditions, we can apply the inductive hypothesis to the subsets and {λ i } n i=1 to conclude that for all i ∈ [n], i = C λi . Then for any i, we have n+1 Example 1. In this example we remark that Proposition 1 can fail to hold if the orthonormal sets are not assumed to be generic: where e 1 , e 2 ∈ R 3 are the first and second standard basis vectors, and let ν 1 , ν 2 , µ 1 ∈ R. Completing x 1 i1 and x 2 i2 to orthonormal bases, we have of this matrix, we construct L L L = 0 0 , which has nullity 2 > 1.
Note that while these are valid counterexamples, L L L having nullity greater than 1 is not an issue in practice since any vector in the nullspace can be used for the symmetrizing procedure.
Thus, we let i be the entries of a non-zero vector in nullsp(L L L ). It follows that the top-left r X1 × r X1 corner blocks of LX 1 TX 1 and Q MQ are equal, where M = CM for some constant C ∈ R. In fact, since multiplication on the left by a diagonal matrix corresponds to scaling of rows, this means that the first r X1 rows of LX 1 TX 1 and Q MQ are equal. Note that the entries in rows r X1 + 1 to n are all 0. By the symmetry of Q MQ , we can find the first r X1 columns of Q MQ . Hence, we have a symmetric S ∈ R n×n with a bottom-right (n − r X1 ) × (n − r X1 ) corner block of 0's The non-zero entries of S and Q MQ are equal. The next step is to determine what entries s ij should be filled in block of 0's in S such that the resulting matrix S is equal to Q MQ . Now we use the DRC: Suppose r X1 ≥ r X2 . Then after performing Gaussian elimination on S, s r X 1 +1,i = 0 for r X1 + 1 ≤ i ≤ n. Thus we replace these entries in S with variables x i and perform Gaussian elimination. This will result in a linear expression for each x i . Setting them equal to 0, we solve for the unique values s r X 1 +1,i . In addition, since S is symmetric, we know the values s i,r X 1 +1 .
We then repeat this procedure n − r X1 more times. The final iteration will be S (n−(r X 1 +1)) = S. Since Q MQ = X 1 X 2 MX 2 X 1 is such a matrix satisfying this Gaussian elimination property when r X1 ≥ r X2 , it follows from the uniqueness ofs ij that S = X 1 X 2 MX 2 X 1 . Thus we can obtain x 2 i2 UTPS via an eigendecomposition on X 1 SX 1 = X 2 MX 2 . If, however, r X2 > r X1 , then we can still apply this procedure to obtain an orthonormal set. In either case, we denote the set as {x 2,LR i2 } r X 2 ,LR i2=1 . This completes the description of the kernel completion method.
may differ from the set {x 2 i2 }, even UTPS. This could happen if r X2,LR > r X1 , since then r X2 ≤ r X3 ≤ ... ≤ r X L is true in the DRC. In this case, the direction of the decomposition (from left to right) is incorrect and we terminate the decomposition. All other sets {x j,LR is the correct set and we continue to sequentially apply kernel completion. For our example, the next iteration of kernel completion to find the set {x 3 i3 } would be applied to where v ∈ R n is a generic vector. We check on each iteration whether r Xj ,LR > r Xj−1,LR and take the necessary course of action. The decomposition finishes by either terminating at some point or obtaining the orthonormal set {x L−1,LR Step 2: Decomposing the train from right to left This step is the same as Step 1 but with "LR" interchanged with "RL" and X 1 interchanged with X L , X 2 interchanged with X L−1 , etc.
Step 3: Choosing the correct orthonormal sets from the two decompositions have different ranks, then whichever has the higher rank is the correct set UTPS. The potential issue is choosing a correct set when r Xj ,LR = r Xj ,RL , for one of the sets might be incorrect. This, however, cannot happen; if, without loss of generality, is incorrect, then it would have been because r Xj > r X j , for some 1 ≤ j < j. But then by the DRC, we would have observed that r Xj = r Xj ,RL > r Xj ,LR , contradicting r Xj ,LR = r Xj ,RL .
Step 4: Finding the coefficients Since the tensor R defined in (3) is a rank-1 tensor, we can apply a tensor decomposition algorithm such as the Alternating Least Squares (ALS) algorithm presented in [3] to find the coefficients. We have observed in simulations that the rank-1 ALS algorithm always converges in one step. The coefficients are, of course, not unique since if ν 1 , ..., ν L ∈ R \ {0} are non-units, then This completes the description of our algorithm for decomposing tensors according to an orthogonal tensor train of length L ≥ 3 with symmetric and orthogonal carriages. The pseudocode can be found in Algorithm 2, and our simulation results are in Section 7. Notice that if L = 1, then Algorithm 2 is none other than the slice method of Kolda [4], and if L = 2, then Algorithm 2 reduces to Algorithm 1 for symmetric orthogonal tensor trains. We have distinguished these two algorithms since the whitening procedure can be applied to symmetric but non-orthogonal tensor trains of length 2. Extending this procedure to longer trains is one of our open problems in Section 8. From our analysis in this section, we have deduced the following result: where λ j ij ∈ R are generic coefficients and {x j ij } r X j ij =1 ⊂ R n are generic orthonormal sets of rank r Xj for all j ∈ [L], and T satisfies the Decreasing Ranks Condition (c.f. Definition 1), then Algorithm 2 recovers all of the ranks of T , and all of the vectors and coefficients of T UTPS.

Algorithm 2 Orthogonal Symmetric Decomposition for Tensor Trains of Length L ≥ 3
Input: T ∈ R n L+2 , an orthogonal symmetric tensor train of length L ≥ 3 satisfying the DRC Output: λ j ∈ R r X j , {x j ij } r X j ij =1 ⊂ R n , the coefficients and orthonormal vectors which generate T , correct UTPS in the vectors and non-zero scaling in the coefficients, for j ∈ [L] A ← n × n orthogonal matrix whose first r X j−1,LR columns are ← eigenpairs of ASA with σ ij = 0 22: if r Xj ,LR > r Xj−1,LR then 23: end if 26: end for 27: Repeat steps 1 to 26, with "LR" interchanged with "RL", j = 2, ..., L − 1 with j = L − 1, ..., 2, j − 1 with j + 1, step 1 replaced with T (v, ..., v, ·, ·), and step 7 replaced with S ← A T A ij =1 , whichever of r X j,LR or r X j,RL is greater, 2 ≤ j ≤ L − 1. If they are the same, make an arbitrary choice 29: R ← R r X 1 ×...×r X L tensor with entries Rî

Orthogonal Decomposition of Tensor Trains of Length 2
We now turn to the problem of decomposing 4-tensors according to tensor trains of length 2 such that the carriages are orthogonal but not necessarily symmetric.
Problem 3. Let T ∈ R n A ×n B ×n D ×n E be a 4-tensor admitting the following decomposition: ⊂ R n F are generic orthonormal sets, λ i , µ j ∈ R are generic, and d = n C = n F (necessarily d ≥ r L , r R due to orthonormality). In other words, assume that T is a tensor train generated by Given T , we wish to find all ranks, vectors, and coefficients.
We first note that the vectors a i , b i , d j , and e j can be found by adapting Kolda's slice method for orthogonal tensors [4]. As before, define the matrices where Σ ∈ R r L ×r L , Γ ∈ R r R ×r R are diagonal with non-zero diagonal entries, and A ∈ R n A ×r L , B ∈ R n B ×r L , C ∈ R n C ×r L D ∈ R n D ×r R , E ∈ R n E ×r R , F ∈ R n F ×r R are the matrices that have {a i }, {b i }, {c i } {d j }, {e j }, and {f j } as their columns, respectively. Then A, B, D, and E can be found by finding the SVD of S L and S R . Next we see that for each i ∈ [r L ], j ∈ [r R ], T (a i , b i , d j , e j ) = λ i µ j c i , f j , which defines the entries of a matrix X ∈ R r L ×r R . Note that n C = n F since they correspond to the same edge of the tensor network. Let d = n C = n F , and suppose we know d. Then we can define the 0-padded d × d matrix where C, F ∈ R d×d are orthogonal matrices whose first r L and r R columns are C and F. To solve Problem 3, we need to find Λ, M, C, and F such that X = Λ C FM. If a solution exists, it is not unique since C = C F and F = I d×d is also a solution. Hence we can instead ask to find an orthogonal matrix Q = C F ∈ R d×d such that X = ΛQM. Lastly, we can assume Λ and M have non-negative diagonal entries since if they do not, we can "push" the negative signs of these entries into Q without affecting its orthogonality. Note that even with this assumption, if a solution exists, then Λ and M are still not unique due to a global non-zero scaling.

Algorithm 3 Orthogonal Decomposition for Tensor Trains of Length 2
Input: T ∈ R n A ×n B ×n D ×n E , a length 2 orthogonal tensor train Output: r L and r R , the ranks of L and R  15: c i ← the first r L rows of Q 16: f j ← the first r R standard basis vectors of R d A matrix of the form ΛQM is said to have a Diagonal-Orthogonal-Diagonal Decomposition (DODD), which is the subject of Section 6. Note that there could be many possible d for which X in (4) has a DODD. Thus to solve Problem 3 (see pseudocode in Algorithm 3), we solve the problem of finding such a d and Λ, Q, and M; this is formally stated in Problem 4. When r L = r R and we know that d = r L = r R gives an X with a DODD, which we call the square case, we show that Problem 4 can be solved using two different approaches: one based on Sinkhorn's algorithm [14] and the other based on the Tandem Procrustes algorithm [15]. We call the complementary case the general case and we present a solution by showing that the Procrustes-based algorithm can be generalized. Lastly, similar to Remark 3.3 in Section 3, we add that one can easily generalize to the case where each carriage is an m-tensor, for some m > 3, but cannot easily address the case where there is more than one contracted edge between L and R.

Applications
The decomposition of length-2 tensor trains is closely related to the tensor hypercontraction (THC) method of compressing the electron repulsion integral (ERI) tensor common in electronic structure theory, first described in [17] and [18]. One seeks an approximate factorization of the fourth-order ERI tensor R of the form where x P are vectors and Z is known as the intermediate matrix. To find this, the PARAFAC-THC method is proposed, which first finds a factorization of R as the contraction of two third-order tensors S and T , and three matrices L, M, and N, as shown in Figure 6. Using the Alternating Least Squares algorithm [3], a factorization is obtained, and setting Z P,Q = A,B,C,D (y P ) A L A,B M B,C N C,D (y Q ) D completes the factorization (5). Now suppose {x P } P is an orthonormal set. Then by fully contracting the tensor network in Figure 6 into a length-2 tensor train, we see that the vectors x P can also be approximated (or in the square case, found exactly) using Algorithm 3. Following this, the values of Z P,Q can easily be recovered without knowing L, M, or N. [19] also proposes another approximate solution using Fast Fourier Transforms. Figure 6: A PARAFAC-THC decomposition [18] of the ERI tensor.

Matrix Diagonal-Orthogonal-Diagonal Decomposition
We now discuss the problem of finding a diagonal-orthogonal-diagonal decomposition of a given matrix. Let m = r L and n = r R , from the previous section. where Q ∈ R d×d is an orthogonal matrix, and λ i , µ j ∈ R. Given such an X, find Λ, Q, and M.

The Square Case: A Sinkhorn-Based Algorithm
We first solve the case when m = n and we know that X has a DODD for d = m = n. We do this by modifying Sinkhorn's algorithm [14], and we additionally require that X has no entries equal to 0. Let · * 2 and * √ · denote the entry-wise square and square root of a matrix, and let denote the Hadamard product of two matrices. Then, X * 2 = Λ 2 Q * 2 M 2 . Every entry of X * 2 is positive, and since Q is orthogonal, Q * 2 is a doubly stochastic matrix. Sinkhorn's algorithm [14] applied to X * 2 first rescales all rows of the matrix X * 2 so they each sum to 1, then it rescales all columns so they sum to 1, and then the rows, alternating until convergence. The total rescaling of rows and columns is recorded in two diagonal matrices, Λ and M, respectively. Sinkhorn's Theorem [14] guarantees that this algorithm converges and yields the unique doubly stochastic matrix Q * 2 and two positive diagonal matrices Λ and M, unique up to non-zero, positive scaling, such that X * 2 = ΛQ * 2 M. Hence, we obtain Q * 2 = Λ −1 X * 2 M −1 . The absolute values of the entries of Q are therefore * Q * 2 . Since Λ and M are positive diagonal matrices, the sign modifications to the entries of * Q * 2 required for it to be orthogonal can be obtained directly from the sign pattern of X. That is, if S ∈ R d×d is the matrix such that S ij = 1 if X ij > 0 and S ij = −1 if X ij < 0, then S * Q * 2 is an orthogonal matrix. Lastly, taking the entry-wise square root of both sides of ΛX * 2 M = Q * 2 together with the sign matrix S, . We present pseudocode in Algorithm 4. From this discussion, we have also obtained the following result: Theorem 3. If X ∈ R n×n admits a decomposition of the form X = ΛQM where Λ, M ∈ R n×n are diagonal matrices and Q ∈ R n×n is an orthogonal matrix, and X has no entries equal to 0, then Algorithm 4 is guaranteed to find Λ, Q, and M, UTPS.

The Square Case: A Procrustes-Based Algorithm
The Tandem Procrustes Algorithm [15] generalizes the well-known Procrustes Problem [20] to orthogonal, but not necessarily orthonormal columns: Given a target matrix A ∈ R m×n and a starting matrix Algorithm 4 Sinkhorn-Based Algorithm for Square DODD Input: X ∈ R d×d , a matrix with an existing DODD Output: Diagonal matrices Λ, M ∈ R d×d and orthogonal matrix Q ∈ R d×d such that X = ΛQM Algorithm: 1: procedure square_dodd(X)

2:
Perform Sinkhorn's algorithm on X * 2 , producing diagonal matricesΛ andM, and a doubly stochastic matrixQ such thatΛX * 2M =Q 3: return Λ, Q, M 9: end procedure B ∈ R p×n , it finds a matrix U ∈ R m×p with orthogonal columns such that ||A − UB|| 2 is locally minimized, where U = VD with diagonal matrix D ∈ R p×p and matrix with orthonormal columns V ∈ R m×p . We recall the Tandem Procrustes algorithm in Algorithm 5.

Input:
A ∈ R m×n , the target matrix B ∈ R p×n , the starting matrix Output: Diagonal matrix D ∈ R p×p and matrix with orthonormal columns V ∈ R m×p such that ||A − VDB|| 2 is locally minimized Algorithm:  We present an approximate solution to Problem 4 in the square case. In our solution, we initialize Q = X, Λ = M = I n . We first call Tandem Procrustes on A = Q and B = I n to find a matrix with orthogonal columns VD that approximates the rows of Q: Q ≈ DV . We update Λ and Q by multiplying Λ with D on the right and Q with D −1 on the left, so that now Q ≈ V . Then we call Tandem Procrustes on A = Q and B = I n to find a matrix with orthogonal columns VD that approximates the columns of Q: Q ≈ VD. We update M and Q by multiplying M with D on the left and Q with D −1 on the right, so that now Q ≈ V. These alternating steps are repeated until the final Q found is an orthogonal matrix. Our solution algorithm is presented in Algorithm 6.
Algorithm 6 Procrustes-Based Algorithm for Square DODD Input: X ∈ R d×d , a matrix with an existing DODD Output: Diagonal matrices Λ, M ∈ R d×d and orthogonal matrix Q ∈ R d×d such that X = ΛQM Algorithm:

Remarks On The Square Case DODD Algorithms
Notice that if a DODD exists for a square matrix, then the Sinkhorn-based algorithm recovers the decomposition only because the entry-wise square of an orthogonal matrix is incidentally doubly-stochastic. If one seeks to find a DODD or a DODD approximation for a square matrix that is not known to have a DODD, then the Sinkhorn-based algorithm can fail. On the other hand, the Procrustes-based algorithm finds a DODD approximation that is locally optimal in this case. We will soon see that this property of the Procrustes-based algorithm allows us to approximately solve the general case for Problem 4, which the Sinkhorn-based algorithm cannot. A numerical and runtime comparison between the Sinkhorn-and Procrustes-based solutions for the square case is presented in Section 7.3.

The General Case
Consider the general case: either m = n or we do not know what value of d ≥ m, n gives a DODD for the d × d 0-padding of X. Based on numerical results in Section 7.3, we conjecture that a DODD always exists for X in Problem 4 for some sufficiently large d ≥ m, n. Hence, suppose we already know the value of such a d. Then we present an approximate solution to Problem 4 where we find Λ, Q, and M. The solution involves modifying Algorithm 6. We can initialize λ i = µ j = 1 and the top-left m × n corner block of Q to be X. The issue we face is how to initialize the other entries of Q, which we will call Q, since the corresponding entries in X are all 0, and the final solution Q must be an orthogonal matrix: Let us call steps 6 to 11 in Algorithm 6 the Tandem Procrustes iterations. We address this issue by initializing all entries of Q to be random numbers and then iteratively correcting these entries, replacing them with entries from an orthogonal matrix, until Q is itself orthogonal. At the same time, we apply the Tandem Procrustes iterations to find Λ and M, thus solving Problem 4. Choose a positive integer , which is to be thought of as a "learning rate". First perform the Tandem Procrustes iterations times. Then replace the entries of Q with the corresponding entries V in V, where V is the most recent orthogonal matrix found by the Tandem Procrustes iterations: where λ i and µ j are the most recent non-zero entries of Λ and M. Finally, repeat these two steps until convergence. In Section 7.3, we demonstrate the numerical performance of this procedure. Also note that as in Algorithm 6, one could have performed the Tandem Procrustes iterations until convergence rather than terminating after times. However, our numerical results suggest that choosing to be small and repeating these two steps greatly decreases the runtime of our solution while producing equally effective solutions. Algorithm 7 presents pseudocode of this procedure.
Algorithm 7 An Algorithm for General DODD Input: X ∈ R m×n , a matrix d ≥ m, n such that the 0-padded matrix X ∈ R d×d , with X in the top-left corner, admits a DODD ≥ 1, the learning rate Output: Diagonal matrices Λ, M ∈ R d×d and orthogonal matrix Q ∈ R d×d such that X = ΛQM Algorithm: .., 1, 0, ..., 0), where 1 appears m times

Numerical Results
The results for symmetric trains of length 2 and the DODD were computed using Python 3.7.6 and the TensorLy package [

Symmetric Orthogonal Tensor Trains of Length 2
Here we tested Algorithm 1. The parameters n, r A , r B (cf. Section 3), control the size of Problem 1, and for each choice of these parameters, 100 tests were conducted. In the same approach to testing as Kolda's slice method [4], a tensor T * with an exact decomposition is artificially constructed. Then a new tensor T is formed by corrupting T * with Gaussian noise: where σ is a noise parameter and N is a tensor of the same dimensions as T * whose entries are independently sampled, standard normally distributed numbers. We tested for σ ∈ {0, 10 −2 , 10 −6 }. If whitening is applied, then the algorithm was given 200 iterations to find PSD matrices C A and C B before declaring failure. The relative error between the solution returned by the algorithm T and T was then computed: For a successful solution, we expect the average log 10 relative error when raised to the power of 10 to be on the order of σ. The results are listed in Table 1. We see that with no noise, the algorithm in both the orthogonal and non-orthogonal cases performed well, with PSD matrices found in most tests for whitening. In the orthogonal case when noise is added, the relative error was generally one order of magnitude more than the noise, similar to Kolda's results in [4]. As in Kolda's tests for the slice method, the poorest performance was for noise in the non-orthogonal case. Kolda hypothesizes that this is due to the noisy tensor having a rank higher than the rank found by the algorithm [4].

Symmetric Orthogonal Tensor Trains of Length L ≥ 3
We tested Algorithm 2 on tensor trains satisfying the DRC for varying n, L, and ranks of carriages r Xi , where i ∈ [L] (cf. Section 4). The method of testing was identical to that of Section 7.1. The results are presented in Table 2. Without noise, the algorithm was successful, but even in the presence of small noise, the relative    error was high, likely due to the rank of the tensor increasing, as Kolda hypothesized. The runtimes were also diverse due to dependence on r Xi (since kernel completion will take longer for lower ranks) and stopping criteria whenever the rank increases.

Matrix Diagonal-Orthogonal-Diagonal Decomposition
Next we tested the Sinkhorn-based Algorithm 4 and Procrustes-based Algorithm 6 for the square case of DODD. We constructed 100 square matrices with exact decompositions for sizes n ∈ {3, 10, 25}. For Algorithm 4, convergence was said to be achieved when the sum of squares of differences between the row and column sums of Q * 2 and 1 was less than 10 −28 . For Algorithm 6, the learning rate was set to 1, and convergence was said to be achieved when the most recent call to the Tandem Procrustes algorithm satisfied ||Q−VD|| 2 < 10 −28 . Both algorithms were given a maximum of 1000 iterations before terminating if convergence was not achieved before then. Both algorithms were said to be successful if the relative error of QQ and I n was within 10 −10 . For all values of n and all tests, both algorithms converged and were successful. Figures 8 and 9 show histograms of the number of iterations before convergence for each algorithm. We saw that as n increased, the number of such iterations generally decreased for both algorithms. The Procrustesbased algorithm generally required fewer iterations than the Sinkhorn-based algorithm. However, as shown in Figure 7a, the overall runtime for Algorithm 4 was faster than that of Algorithm 6, since the Tandem Procrustes algorithm requires computing polar factors. Figure 7b shows a sample of the error trajectory of the Tandem Procrustes algorithm on a random 10 × 10 matrix. Even for a matrix unlikely to have a square DODD, the algorithm appears to converge to a local minimum.
Lastly, we examined Algorithm 7 for the general case of DODD. For m = 6, n = 5, and each value of d ≥ m, n tested, we constructed 100 m × n matrices whose 0-padded d × d matrix admitted exact decompositions ("exact"), and 100 m × n matrices whose entries are independently sampled from a normal distribution with mean 0 and variance 5 2 ("non-exact"). We then ran Algorithm 7 on these matrices, inputting d for both exact and non-exact tests, with learning rate = 2. The convergence and success criteria was the same as for the square case. The results are presented in Table 3. We saw that when d = max{m, n}, the algorithm mostly performed well for the exact tests and poorly performed for the non-exact tests. What was intriguing, however, was that for d = {m, n}+1, we observed no successes for both tests, and as d increased, the frequency of success also increased, reaching complete success even for the non-exact tests when d was sufficiently large. This lead us to the following conjecture: Conjecture 1. The d×d 0-padding of any matrix X ∈ R m×n admits a DODD, for some d ≥ m, n sufficiently large. This means that every matrix X ∈ R m×n has a decomposition of the form X = ΛQM, where Λ ∈ R m×m and M ∈ R n×n are diagonal matrices and Q ∈ R m×n is the m×n top-left corner block of an orthogonal matrix Q ∈ R d×d . Equivalently, X has a decomposition of the form X = ΛC FM, where C ∈ R m×d is a matrix with orthonormal rows and F ∈ R d×n is a matrix with orthonormal columns (we can choose C to be the first m rows of Q and F to have the first n standard basis vectors in R d as columns).  Table 3: General case DODD for both the "exact" and "non-exact" case. For each permutation of control parameters, 100 tests are run, with each test given a maximum of 1000 iterations to converge.

Conclusion and Future Work
In this paper, we studied the decomposition of tensors into tensor trains whose carriages are symmetric or orthogonally decomposable tensors. For length-2 trains, we showed that we can find the decomposition when the carriages are symmetric and odeco (Section 3.1), only symmetric (Section 3.2), or only odeco (Section 5). For longer length tensor trains, we studied the case when the carriages are both symmetric and odeco. For such networks, we provided algorithms for decomposition which require the Decreasing Ranks Condition (Section 4). In Section 7, we provided numerical results that support our findings.
A variety of open problems arose during our study. We showed that Kolda's whitening procedure for symmetric tensors could be adapted to the similar case of a tensor train of length 2, whereby a tensor with a symmetric decomposition could be transformed into an orthogonal problem. It remains an open problem to find an efficient equivalent to longer tensor trains, and so encompass a greater family of tensors for such a decomposition.
Problem 5. Can we find a whitening procedure for symmetrically decomposable tensor trains of length greater than 2?
One of the interesting mathematical problems that arose during our study was that of finding a Diagonal-Orthogonal-Diagonal Decomposition for matrices (Section 6). We are not aware of any linear algebra literature that has previously studied this topic, and there are a number of issues to still address.

Problem 6.
For what hidden dimension d can we decompose any given matrix via a Diagonal-Orthogonal-Diagonal Decomposition?
We saw via our numerical tests that taking d to be the maximum of the two dimensions of a matrix and applying our algorithm did not result in an exact DODD, but rather convergence to a suboptimal approximation. When d, however, was strictly greater than both dimensions, we saw that the error output at each iteration of our algorithm continuously decreased, and given either a sufficiently large maximum threshold of iterations, or given a sufficiently large d for efficient convergence, an exact solution was generally achieved. It remains an open problem to find exactly the minimum such hidden dimension d relative to the dimensions of a given matrix such that an exact DODD exists.
Assuming an exact solution exists, we solved the square case for the DODD, i.e., the case where the hidden dimension d is equal to both dimensions of the given matrix. The question remains for exactly what family of matrices this is the case. Problem 7. Describe the set of matrices (as a set cut out by polynomial equations) that have a square Diagonal-Orthogonal-Diagonal Decomposition.
More broadly, we began our study with a focus on the decomposition of orthogonal tensor networks in general. As a specific network form with many modern applications, and an ability to describe large families of tensors efficiently, tensor trains were the primary case we considered, and due to fruitful progress resulted in being the main focus of this paper. However, we remain interested in other, more general tensor networks.
Problem 8. Can we decompose orthogonal tensor networks for any general tensor network? For example, can we find an orthogonal tensor ring decomposition [16]? Can our methods for orthogonal tensor trains be extended to orthogonal tensor trees?
The primary motivation for this work is to find a structured decomposition that applies to all tensors. Orthogonal tensor train decompositions do not apply to all tensors of a given size. As discussed earlier, the rank of a tensor decomposing this way can be up to O(n d−2 ) while the generic rank of a dimension-n order-d tensor is O(n d−1 ). Orthogonality will, of course, make the orthogonal tensor train decomposition even more restrictive. So, a natural question is which orthogonal tensor networks can represent all tensors of a given size. Another direction for future work is to devise a method which approximately decomposes a given tensor as a tensor train.
Problem 9. What are the simplest orthogonal tensor networks according to which any tensor of a given size can be decomposed? Given a tensor and a tensor network diagram, can we approximately decompose the tensor according to that diagram?
Along the lines of the previous question, it would be interesting to give an implicit description of the set of tensors that decompose according to a given orthogonal tensor network.
Problem 10. Find the polynomial equations that define the set of tensors that decompose according to a given orthogonal tensor network. Is the set of such tensors Zariski closed?
Finding the eigenvectors and singular vector tuples of a general tensor is an NP-hard problem [2]. However, the problem is easy for the family of orthogonally decomposable tensors [7,8]. How about for orthogonal tensor networks? Problem 11. Given a tensor decomposing according to an orthogonal tensor network, can we find its singular vector tuples efficiently? Can we describe the set of all of its singular vector tuples in terms of its orthogonal tensor network decomposition?
Answering these questions would provide fundamental progress towards finding a structured tensor decomposition for any tensor, analogous to the singular value decomposition for matrices.