1 Introduction↩︎

Tensors are multidimensional extensions of vectors and matrices [1], used pervasively both for modelling physics (continuum mechanics, electromagnetism, general relativity, quantum mechanics, ...), and as a tool for computer science (computer vision, linear system solving, machine learning, quantum computing, ...).

A tensor is a rather big object, the amount of data required to store it depending exponentially in its order (i.e. the number of dimensions; matrices and vectors being respectively order-2 and order-1 tensors). To deal with this issue, schemes have been devised to compress the data, by decomposing the tensor into a tensor network: a collection of smaller-ordered tensors composed together following a given topology [2]. It is then possible to exploit the Singular Value Decomposition (SVD) to decrease the dimensions between the tensors, either exactly (by exploiting the rank-revealing property of the SVD), or approximately, by chopping off singular values below a given threshold, and hence performing low-rank approximations [3].

A first caveat of the SVD is that it is in general not unique, which means that SVD-based reduction schemes never enjoy a full uniqueness property [1]. The second is that it is only applicable to matrices over \(\mathbb{R}\) or \(\mathbb{C}\) (hence using floats in practice). For most considered applications, like the aforementioned, this is fine. However, it recently was shown that tensor trains (tensor networks consisting of order-3 tensors connected in a line) over \(\mathbb{F}_2\) could be used to generalise Binary Decision Diagrams (BDDs) [4], [5]. In particular, it has been shown that there exist families of tensor trains over \(\mathbb{F}_2\) of polynomial size, that would require BDDs of exponential size to represent them [4].

Binary decision diagrams [6] are a data structure used to represent boolean functions, and that make use of the functions’ structure to reduce their size. These diagrams and their numerous variants (Zero-suppressed DDs [7], BDDs with complemented edges [8], ...) can be and have been used for a plethora of applications, such as model checking [9], model counting [10], circuit synthesis [11], regular expression matching [12], multiobjective discrete optimization [13], ..., or simply for the storing of boolean polynomials [14].

The fact that BDDs are special cases of tensor trains over \(\mathbb{F}_2\) gives an interesting bridge between two rather distant fields. In particular, it opens the possibility to use tensor trains rather than BDDs for many of the latter’s application, with potentially better outcomes due to the better compactness of tensor trains. However, to properly assess the advantage obtained by using tensor trains, one would first need to address the two aforementioned caveats of SVD-based tensor train reductions: application in \(\mathbb{F}_2\), and uniqueness.

We address this in this paper, by providing a unique normal form for tensor trains over arbitrary fields, that minimises their size, and can be reached in polynomial time in their original size. We also provide an algorithm that directly builds a tensor train in normal form out of any arbitrary full tensor.

The developments in this paper are based on a matrix decomposition that is different from the SVD: the LU decomposition, whereby a matrix is decomposed as the product of a lower triangular matrix and an upper triangular matrix. Although a matrix may not have an LU decomposition, this can be fixed by first permuting its rows, leading to the so-called PLU decomposition. While this one always exists, it is not unique in general, which is inconvenient for our purposes. We hence use a variation of the LU decomposition, the LDPU decomposition, that we show to be unique. Crucially, just like the SVD, the LDPU decomposition is rank-revealing, the central feature that allows us to reduce the size of the tensor trains. While the LDPU (and subsequent tensor train normal form and reduction strategy) theoretically work with any field, in practice, it may suffer from exponential bit complexity and numerical instability, hence it is best fitted for finite fields, where such concerns disappear.

The paper is structured as follows. In 2, we give some necessary background on tensor trains, a short explanation of the link between these and decision diagrams, and a string-diagram based representation to simplify their manipulation. In 3, we define the LDPU decomposition, show its uniqueness, and give a few useful lemmas about matrices in row-echelon form. In 4, we define the normal form, as well as an algorithm that directly turns a full tensor into its normal form as a tensor train. The existence and uniqueness of the tensor train is then shown there. Finally, in 5, we define the poly-time reduction strategy that brings any tensor train into the unique normal form defined previously. We also show how the first part of the reduction is already enough to infer some useful properties about the tensor train, namely, whether it is null or not (which induces an algorithm to check equivalence between two tensor trains), or what are its leading index and leading coefficient. All proofs are given in full in the appendix.

2 Tensor Trains↩︎

Throughout this paper, we will work with tensors over arbitrary fields. Let’s fix one, that we denote \(F\). We write \(0\) and \(1\) for its additive and multiplicative neutral elements. For any tensor \(T\in F^{n_1\times ... \times n_d}\), we write \(T[i_1]\) for the tensor \(F^{n_2\times ... \times n_d}\) obtained by setting the first index of \(T\) to \(i_1\). Hence, for a matrix \(A\), \(A[i]\) corresponds to the row \(i\) of \(A\), and for a vector \(v\), \(v[i]\) corresponds to the \(i\)-th element of \(v\). We write \(T[i_1,...,i_k]\) for \(T[i_1][i_2]...[i_k]\).

The notion of tensor trains is sometimes generalised to the case where \(r_0=r_d\) can take any value, which generalises the interpretation and gives more leeway in the manipulation. However, the more constrained version \(r_0=r_d=1\) makes the reduction of tensor trains much easier, in particular because reaching the minimum in TT-ranks can be made one TT-rank at a time (the TT-ranks reach a minimum iff no individual TT-rank can be reduced).

To make the presentation simpler in the following, we allow modes of tensors to be of dimension \(0\). This also means we allow matrices with \(0\) rows or with \(0\) columns. Since such tensors have no entries, they are, by convention, considered null. Tensor trains allow in certain cases a more compact representation of the full tensor. That representation is polynomial in \(d\), \(n\) and \(r\), and so are many operations on the tensor trains [15]. For instance, building a tensor train that represents a linear combination of tensor trains, the element-wise product (a.k.a Schur product) or the Kronecker product of tensor trains, can be done in polynomial time (and space). Many operations that we may want to apply on a tensor may also be implemented as a poly-time operation on the tensor train. There are hence many occurrences where the full tensor is never computed: the computation starts with a standard tensor of which we know an efficient tensor train representation, and then all operations are directly applied on the tensor train. It is then paramount to have efficient reduction strategies to keep the size of the tensor train under control.

2.1 BDDs as special cases↩︎

As a main motivation for the development of tensor trains for arbitrary fields, we briefly explain how Binary Decision Diagrams (BDDs) are special cases of tensor trains over \(\mathbb{F}_2\), the \(2\)-element Galois field; where the cores are \(1\)-row-sparse [4]. A BDD (without variable skipping¹) is a rooted directed acyclic graph, where the nodes can be partitioned into \(h\) groups such that edges can only go from one group to the next. The last group consists of 2 nodes labelled \(\bot\) and \(\top\). For each node except the last two, there are exactly two outgoing edges, labelled \(0\) and \(1\) respectively. Each group of nodes (except the last) is assigned a distinct variable, which allows the structure to represent a boolean function over these variables: one can easily evaluate the function for a given variable assignment, by following the path where the edge labels correspond to the values assigned to the variables. The evaluation of the function is then \(1\) if and only if the path ends up in the node labelled \(\top\).

Such a BDD is easy to translate into a tensor train \(T\in \mathbb{F}_2^{2\times...\times 2}\). Indeed, consider the (ordered) groups \(i\) and \(i+1\). These nodes together with the \(0\)-edges between them form a bipartite graph, which we can represent by an \(\mathbb{F}_2\) biadjacency matrix (with the rows corresponding to the nodes of group \(i\), and the columns those of group \(i+1\)). For each group \(i\) (but the last), we can then easily build a core \(C_i\) such that \(C_i[0]\) is the biadjacency matrix for the \(0\)-edges from group \(i\) to group \(i+1\), and \(C_i[1]\) is the biadjacency matrix of the \(1\)-edges from group \(i\) to group \(i+1\). For the penultimate group, we only look at the edges to the \(\top\) node (i.e. we treat the last group as if it only contained the \(\top\) node). It has been proven that this construction yields a tensor train that encodes the same boolean function. An example of a tensor train derived from a BDD is given in 1.

It is clear now that all matrices in the cores obtained directly from a BDD are \(1\)-row-sparse, for each node has at most \(1\) outgoing \(0\)-edge, and at most \(1\) outgoing \(1\)-edge. While TTs can be exponentially more succinct than BDDs, BDDs can at best be polynomially more succinct than TTs, thanks to the BDD to TT translation that we just described [4].

By extension, Weighted Binary Decision Diagrams (WBDDs) [16], [17] are sparse special cases of tensor trains over the corresponding field. Notice that just like for variable ordering in BDDs [18], [19], core ordering in tensor trains may dramatically change the size of the data structure [1], [20].

Since arbitrary matrices can be seen as biadjacency matrices for associated bipartite graph, we can still have a DAG-based graphical interpretation of arbitrary tensor trains. In this case, every non-terminal node can have arbitrarily many outgoing \(0\)-edges et \(1\)-edges, and each node can be seen as encoding a vector (or tensor), built by concatenation of linear combinations of the vectors encoded by the nodes in the next group (the children nodes), starting at the \(\top\) node interpreted as \(\begin{pmatrix}1\end{pmatrix}\). While in BDDs (resp. WBDDs), nodes can be merged when they represent equal (resp. colinear) vectors, the reduction defined below can be interpreted as: whenever a node represents a vector that is a linear combination of the others, it can be removed (and connections split amongst the others).

2.2 Graphical representation↩︎

Tensor trains, also called Matrix Product States (MPS) [21], are a special kind of a more general data structure called tensor network [22], [23]. As such, they can be conveniently represented graphically. To keep things completely formal, we will work with string diagrams representing the strict compact-close symmetric monoidal category \(\mathbf{Mat}_F\) whose objects are tuples of natural numbers and morphisms \(T:F^{n_1}\otimes...\otimes F^{n_k}\to F^{m_1}\otimes...\otimes F^{m_\ell}\) are linear maps from \(F^{n_1}\otimes...\otimes F^{n_k}\) to \(F^{m_1}\otimes...\otimes F^{m_\ell}\) [24], [25].

The tensor product of two vector spaces \(F^n\otimes F^m\) is another vector space, that we can define as follows: Each space \(F^{n}\) has a canonical orthonormal basis \(\mathcal{B}_n = \{e_i\}_{0\leq i < n}\) where \(e_i\) is the length-\(n\) vector with entry \(1\) at index \(i\) and \(0\) everywhere else. In that case, \(F^n\otimes F^m\) is the vector space whose canonical basis is \(\{e_i\otimes e_j\}_{\substack{0\leq i< n\\0\leq j < m}}\). This generalises to multiple tensors of the base vector spaces \(F^n\). Since the category is considered strict, we consider the associativity of \((-\otimes-)\) as an equality. If \(T:F^{n_1}\otimes...\otimes F^{n_k}\to F^{m_1}\otimes...\otimes F^{m_\ell}\) and \(S:F^{p_1}\otimes...\otimes F^{p_i}\to F^{q_1}\otimes...\otimes F^{q_j}\) are linear maps, \(T\otimes S\) is a linear map, defined on the basis elements as: \[\begin{align} T\otimes S : F^{n_1}\otimes...\otimes F^{n_k}\otimes F^{p_1}\otimes...\otimes F^{p_i} & \to F^{m_1}\otimes...\otimes F^{m_\ell}\otimes F^{q_1}\otimes...\otimes F^{q_j}\\ e_{r_1}\otimes ... \otimes e_{r_k}\otimes e_{s_1}\otimes ... \otimes e_{s_i} &\mapsto T(e_{r_1}\otimes ... \otimes e_{r_k})\otimes S(e_{s_1}\otimes ... \otimes e_{s_i}) \end{align}\] The tensor product interacts with the categorical composition \((-\circ-)\) according to the bifunctorial law, as follows: if \(S_i\) and \(T_i\) can be composed as \(S_i\circ T_i\), for \(i\in\{0,1\}\), then \((S_0\otimes S_1)\circ(T_0\otimes T_1) = (S_0\circ T_0)\otimes (S_1\circ T_1)\).

There is an isomorphism \(\iota_{n,m}:F^n\otimes F^m\to F^{nm}\), defined on the canonical basis as \(e_i\otimes e_j\mapsto e_{mi+j}\). The isomorphism \(\iota_{n,m}\) is crucial to tensor manipulation, it allows us to reshape tensors to our will. One can easily check that \(\iota_{nm,p}\circ(\iota_{n,m}\otimes I_p) = \iota_{n,mp}\circ(I_n\otimes\iota_{m,p})\), which allows us to unambiguously define a generalised isomorphism \(\iota_{n_1,...,n_k}:F^{n_1}\otimes...\otimes F^{n_k}\to F^{n_1...n_k}\).

The space \(F^1=F\) deserves some additional attention: it is the tensor unit (\(F\otimes F^n\cong F^n\cong F^n\otimes F\)), a property captured by \(\iota_{1,n}\) and \(\iota_{n,1}\). This allows us to see a size-\(n\) column vector as an \(F\to F^n\) morphism, and a size-\(n\) row vector as an \(F^n\to F\) morphism.

The other special morphisms in the category are the ones that make it compact-close and symmetric:

• the symmetry \(\begin{aligned}[t] \sigma_{n,m}: F^n\otimes F^m &\to F^m\otimes F^n\\ e_i\otimes e_j &\mapsto e_j\otimes e_i \end{aligned}\)

• the compact structure \(\begin{aligned}[t] \eta_{n}: F &\to F^n\otimes F^n\\ x &\mapsto x\sum_{i=0}^{n-1}e_i\otimes e_i \end{aligned}\quad \text{and}\quad \begin{aligned}[t] \epsilon_{n}: F^n\otimes F^n &\to F\\ e_i\otimes e_j &\mapsto \left\langle e_i | e_j \right\rangle \end{aligned}\)

This category comes with a few additional axioms, that, for simplicity, we will express in graphical form in the following.

This category is represented with string diagrams using the following notations:

• Identity \(I_{n}:F^n\to F^n\):\(\tikzfigc{14}\)

• Arbitrary tensor \(T : F^{n_1}\otimes ... \otimes F^{n_k} \to F^{m_1}\otimes ... \otimes F^{m_\ell}\):\(\tikzfigc{03}\)

• Symmetry \(\sigma_{n,m}:F^n\otimes F^m \to F^m\otimes F^n\):\(\tikzfigc{07}\)

• Compact structure: \(\eta_n:F\to F^n\otimes F^n:\quad\tikzfigc{05}\)and \(\epsilon_{n}:F^n\otimes F^n\to F:\quad\tikzfigc{06}\)

• Isomorphism \(\iota_{n_1,...,n_k}:F^{n_1}\otimes ... \otimes F^{n_k}\to F^{n_1...n_k}\):\(\tikzfigc{04}\)

• Compositions:\(\tikzfigc{10} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{11}\) \(\tikzfigc{12} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{13}\)

In the following (and already in the compositions above), we way drop the wire annotations to avoid clutter. In the special case of \(I_1\) the identity on \(F\), we may represent it with a dashed line. As particular cases of tensors, a matrix \(M\), a column vector \(v\) and a row vector \(v^\intercal\) are represented respectively as follows: \(\tikzfigc{00}\), \(\tikzfigc{01}\) and \(\tikzfigc{02}\). The compact structure allows us to define the transpose of tensors, as: \(\tikzfigc{15} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{16}\)
Interestingly, the inverse of the isomorphism \(\iota\) is its transpose. We hence may define it as follows: \(\tikzfigc{08}~~:=~~\tikzfigc{09}\). This allows us to define the row vectorised form of an arbitrary tensor:

. This also allows us to define the Kronecker product of two matrices \(M\) and \(N\) can as a reshaping of their tensor product: \(\tikzfigc{00}\)

The graphical notation makes the bifunctorial law a tautology, as each side of the equation is represented in the exact same way as string diagrams. The additional axioms of compact-close symmetric monoidal categories (satisfied by \(\mathbf{Mat}_F\)) are the following, with \(T\) an arbitrary tensor: \[\begin{align} &\qquad\qquad \qquad\qquad \\[1em] &\qquad\qquad \end{align}\] A crucial result on string diagrams for compact-close symmetric monoidal categories, is that one can deform any diagram at will, using the above axioms:

These axioms don’t account for the isomorphism \(\iota\), which is specific to the category \(\mathbf{Mat}_F\). We hence introduce the following additional equations to axiomatise it: \[\begin{align} \qquad\qquad \end{align}\] It is fairly easy to check that these are indeed satisfied. We do not claim any form of completeness here, and simply presented the equations that are necessary for the following.

Within this framework, we may redefine what a tensor train decomposition is. To simplify the presentation, we will work on states. Notice that thanks to the compact closure, tensor train decompositions for states naturally generalise to tensor train decompositions for arbitrary tensors.

Notice that \(C_1\) and \(C_d\) are represented differently, because \(r_0=r_d=1\). We could have represented them as \(None\) and \(None\) respectively to keep them consistent with the rest, albeit at the cost of a more cluttered form.

Before moving on, we want to emphasise the importance of this graphical framework for the following: it is completely formal, and it makes the handling of tensors, and all proofs that manipulate them, much easier and cleaner to present. Without them, proofs can quickly look like a mess with large amounts of indices to properly keep track of during reshapes.

2.3 Reduction↩︎

Tensor train decompositions are clearly not unique, since, for any matrix \(A\) with a left inverse \(A^{-1}\), we have: \[\begin{align} \end{align}\] by pushing \(A\) in \(C_2\) and \(A^{-1}\) in \(C_1\).

In order to reduce as much as possible the memory size required to represent a tensor by a tensor train, we want to find a decomposition that minimises the TT-ranks. One way to do this is to exploit the property above, using rank factorisations of the cores. Such a decomposition has been identified as the Hierarchical Singular Value Decomposition (HSVD). Importantly, such decomposition can be carried out on an existing tensor train, with a complexity that is polynomial in \(d\), \(r\) and \(n\).

While the SVD is numerically stable, provides the best low-rank approximations of arbitrary matrices, and is amenable to several optimisation techniques from numerical methods [27], the HSVD suffers from a few drawbacks:

• It is not unique, although several variations have been introduced to constrain the TT form [21], [28], [29]

• Singular Value Decompositions are only defined on \(\mathbb{R}\) and \(\mathbb{C}\), and can hence only be approximated in practice

• There is no known efficient way of computing the first (or really any) index that contains a value with amplitude \(>\epsilon\) for a given threshold \(\epsilon\geq0\)

The rest of the paper is devoted to providing a new tensor train decomposition that is valid for any field-valued tensor, and reaches optimal TT-ranks in the exact computation case. This decomposition is based on a simpler matrix factorisation than the SVD: the row-echelon form, allowed by mere Gauss-Jordan elimination. As we will see, the tensor train decomposition we obtain is unique, and can be reached from any tensor train in time polynomial in its size.

3 LDPU Decomposition↩︎

In the following, we will use the usual notion of row-echelon form (REF) of a matrix, whereby for every non-zero row, its leading index is (strictly) larger than the leading indices of the non-zero rows above². The row-echelon form process may easily give a rank decomposition of the initial matrix, by removing the zero rows in the echelon form, and the corresponding columns in the transformation matrix that yields the matrix in REF. Since in such a matrix, the rows are linearly independent, we say that such a REF matrix is full row-rank. If the leading coefficients of all the non-zero rows of a matrix in REF are \(1\), we say that it is in unit REF. Obviously, we say that a matrix is in (unit) column-echelon form (CEF) if its transpose is in (unit) row-echelon form.

3.1 Algorithm and uniqueness↩︎

A standard way to compute a REF out of a matrix is to compute a PLU decomposition [30]. This one is however not unique, hence we will use a much less common variation, called LDPU decomposition (which can be partially found in [31] for square matrices). To do so, we need an extra notion on column-echelon form matrices:

Notice that since \(L\) is in column-echelon form, we already have \(i<j\implies L[r_i,j]=0\). Being \(P\)-conditioned hence only adds further constraints when \(i>j\) and \(\sigma^{-1}(j) > \sigma^{-1}(i)\), i.e. when \(\sigma^{-1}\) inverses \(i\) and \(j\).

Notice that since \(r\) is the rank of \(A\), \(U\) has to be full row-rank (i.e. the rows of \(U\) form a linearly independent family of row vectors), and similarly \(L\) has to be full column-rank. Notice also that when working in \(\mathbb{F}_2\), \(D\) is merely the identity matrix, and can hence be ignored.

Obviously, an LDPU decomposition directly provides an LPDU decomposition, by pushing the diagonal matrix through the permutation matrix.

It is fairly easy to get an LDPU decomposition of any matrix, by adapting the standard Gaussian elimination algorithm. First we make sure that all the rows of \(A\) are either null or have unique pivots, by performing as few operations on the matrix as possible: when considering a row, we add another previous row to it only if its leading index is a previously identified pivot. After setting all leading coefficients to \(1\), we get a matrix which is in REF up to a permutation of the rows, which will constitute the permutation matrix \(P\). Since we only add given rows to rows below them, the obtained transformation matrix \(L\) is lower triangular. Moreover, a row is never added to a row with a smaller leading index. This is captured by \(L\) being \(P\)-conditioned. We can then make it unit by extracting the diagonal coefficients, that we store in the diagonal matrix \(D\). Finally, for every zero row in the REF matrix, we remove it, as well as the corresponding row/column in the other three matrices. The obtained algorithm is more formally expressed in [alg:LDPU], where \(\mathbf{0}_{n\times m}\) represents the zero \(n\times m\) matrix.

Notice that if the rank of the matrix is \(0\) (i.e. if the matrix is null), the algorithm still produces a result, with \(L\) having no column, \(U\) having no row, and \(D\) and \(P\) being \(0\times0\) matrices. The existence of this algorithm ensures the existence of an LDPU decomposition for any matrix. This decomposition also turns out to be unique:

As usual, everything is expressed here as a matrix, to simplify the presentation, although in practice, some elements may be more compactly stored (in particular the diagonal and permutation matrices). We are not concerned here with the efficiency and implementation details of the algorithm, for we are more interested in the theoretical results that the decomposition allows.

3.2 Properties of row-echelon forms↩︎

In the following, we will need several results about matrices in row-echelon form, that we state here. The first state that the REF is preserved by Kronecker product:

Notice that the identity matrices are full row-rank unit REF matrices, so the lemma applies as a particular case when either \(U_0\) or \(U_1\) is replaced by the identity. The next lemma states that the REF is also preserved by composition:

The next lemma addresses the existence of right inverses for matrices in REF. It is mostly used to prove the two lemmas afterwards.

The following lemma refines 2 on products of matrices in REF, for the full-rank case:

The final lemma will be necessary for the proof of uniqueness of the normal form, in particular to propagate the uniqueness of individual cores:

4 The Normal Form↩︎

We start by giving an algorithm that builds a tensor train from an arbitrary tensor for which we have the full description, using the LDPU decomposition described above. Since the algorithm is deterministic, we will define its result as the (unique) normal form of the tensor train decomposition. Defining the normal form from this algorithm might seem unorthodox, but we feel like it is simpler to understand the normal form by trying to build a TT directly from the full tensor. We will later show how to reach that normal form from any given tensor train, without the need to compute the full tensor.

4.1 Breaking down arbitrary tensors↩︎

To do so, we need three helper functions that essentially reshape the cores/matrices:

• \(\operatorname{split}(A,n)\) is a “partial vectorisation”, it assumes matrix \(A\) has \(kn\) columns for some \(k\), and breaks each row of \(A\) into \(n\) pieces of size \(k\), to build a block of those \(n\) smaller rows in order, matricially and graphically: \[\renewcommand{\arraystretch}{0.8} \setlength{\arraycolsep}{2pt} \operatorname{split}\left[ \begin{pmatrix} a_{1,1}\cdot \ldots \cdot a_{1,n}\\ \vdots\\ a_{\ell,1}\cdot \ldots \cdot a_{\ell,n} \end{pmatrix} ,n \right] = \scalebox{0.8}{ \begin{pmatrix} a_{1,1}\\ \vdots\\ a_{1,n}\\ \scalebox{1.3}{\boldsymbol{\vdots}}\\ a_{\ell,1}\\ \vdots\\ a_{\ell,n}\\ \end{pmatrix}} \] where \(a_{i,j}\) are row-vectors of size \(k\), and \((\_ \cdot\_)\) is the concatenation of row vectors.

• \(\operatorname{untwine}(L,n,P)\) assumes \(L\) has \(nk\) rows and \(\ell\) columns for some \(k\) and \(\ell\), and that \(P\) is a permutation matrix of size \(k\). The result is a core \(C\in F^{n\times k\times \ell}\) such that \(PC[i] = L[j=i \bmod n]\), where \(L[j=i \bmod n ]\) is matrix \(L\) where we only keep the rows whose indices verify the condition. Graphically: \[None\]

• \(\operatorname{intertwine}\) is the reverse operation: \(\operatorname{intertwine}(\operatorname{untwine}(L,n,P),P)=L\), graphically: \[None\]

The algorithm is given in [alg:NF-from-full-tensor]. For a full tensor \(T\in F^{n_1\times ... \times n_d}\), it is to be called on the row-vectorised version of \(T\) (described as \(\vec{T}^\intercal\) above), the list \([n_1,...,n_d]\), an empty list for the accumulator, and the trivial permutation matrix \(\begin{pmatrix}1\end{pmatrix}\).

We can graphically interpret this algorithm using the tensor network notation, in 2, using \(P_d=U_d=\begin{pmatrix}1\end{pmatrix}\). This graphical understanding will prove useful for the upcoming proofs.

An example of the application of the algorithm can be found in 3, with steps detailed in the appendix.

4.2 Unique characterisation↩︎

The tensor train obtained from this algorithm has two main properties, captured by the following notion of normal form:

Notice that these conditions can be checked in polynomial time in the size of the tensor train.

Notice also that for \(i=1\), the first condition tells us that either \(r_1=0\), or \(C_1\) is non-null, for \(C_1[k]\) is a row vector for all \(k\), hence \(\begin{pmatrix}C_1[0]\mid \ldots\end{pmatrix}\) is obviously in row-echelon form unless it is null. Similarly, for \(i=d\), the second condition is obviously met unless \(C_d\) is null.

Graphically, the two conditions can be expressed as:

• \(None\) is in row-echelon form (unit for \(i>1\)) and full row-rank

• \(None\) is in column-echelon form, full column-rank, and \(P_i\)-conditioned

We call this form normal, for, as we will see, it can be obtained by the upcoming reduction of tensor trains.

It is then possible to show the main result of the paper:

5 Reduction of Tensor Trains↩︎

We show in this section how to turn any tensor train in normal form, by only performing operations on the tensor train itself, not on the full tensor, resulting in an algorithm that is polynomial in the size of the tensor train, and not of the full tensor. This can be seen as the equivalent for tensor trains of the usual (weighted) BDD reductions.

Just like the HSVD form, the algorithm will happen in two sweeps, the first from right to left, and the second from left to right, but the two will be slightly different (although at their core they perform a rank-factorisation to reduce the size of the tensor train). We will show that the sweeps cannot augment the size of the tensor train, and hence that normal form minimises the size of the tensor train (for the given ordering of the modes). The first and second sweeps essentially enforce that the tensor train verifies the first and second normal form conditions, respectively.

5.1 First sweep↩︎

The first sweep is the simplest: for each core starting from the last, we make sure that \(None\) is in unit full row-rank row-echelon form by computing its row-echelon form³: \[\begin{align} \tikzfigc{00} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{01} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{02} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{03} \end{align}\] The new core \(C_i'\) is then a reshape of \(U\), while \(C_{i-1}\) is updated by pushing the matrix \(LDP\) into its input. It should be clear that this operation does not modify the full tensor. The algorithm for the first sweep is given in more detail in [alg:sweep-1].

While this first sweep in itself is not enough to reach the normal form (nor is it enough to get to the minimal size of the tensor train), it is enough to answer the nullity of a tensor train:

This result can be useful for several purposes, without the need to get to the unique normal form. Indeed, one can check if two tensor trains are equal by subtracting one from the other (done in polynomial time), apply a first sweep to the result and check if it is null. When working in \(\mathbb{F}_2\), one can check the satisfiability of a formula by building a tensor train representation of it, and checking for its nullity in the same manner.

Related to satisfiability, one may want to enumerate, or at the very least extract a solution, to the satisfiability of a formula when it is indeed satisfiable. More generally, it is possible to get the index and value of the first non-zero element of a tensor train obtained from the first sweep.

Notice that although the algorithm computes a coefficient of the tensor (the leading one), it does so in \(\mathcal{O}(dnr)\), which is often better than the \(\mathcal{O}(dr^2)\) required to get a coefficient at an arbitrary index – by doing \(d\) matrix-vector products. This is thanks to the echelon structure of the cores, which allows us, for the leading coefficient, to avoid unnecessary computations.

5.2 Second sweep and full reduction↩︎

Although very useful, the first sweep does not yield a unique normal form, nor does it even minimise the TT-ranks.

The second sweep fixes this issue, and is similar in essence, but may seem less intuitive because of the permutations that we have to keep track of, although the graphical notation should make it clearer. We start this time from the first core with trivial permutation \(P_0\) on one element, and enforce the core’s compliance with the second condition of normal forms, by updating the core as follows, assuming permutation \(P\): \[\begin{align} \tikzfigc{00} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{01} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{02} ~~ \underset{\substack{}}{=} ~~ \tikzfigc{03} \end{align}\] The core \(C_{i+1}\) is then updated by pushing \(U\) into its first mode. we then continue with the next core, using the newly computed \(P'\) as the permutation matrix. The algorithm for the second sweep is given in detail in [alg:sweep-2].

Notice that if at some point a null core is created, the nullity is spread from one core to the next during the sweeps. When applying the algorithm to a tensor train representing a zero tensor, after the 1st sweep, the first core will be null, by [prop:nullity-checking]. The second sweep then makes every core null, and the TT-ranks (except \(r_0\) and \(r_d\)) end up being \(0\).

Since the size of a tensor train is smaller that \(dnr^2\) with \(d\) number of modes, \(n\) its largest mode, and \(r\) its largest TT-rank, we can re-express more roughly the runtime of the algorithm as \(\mathcal{O}(dnr^3)\). To compute it, we assumed naive matrix-matrix multiplication, as well as the unoptimised LDPU decomposition described above. If the usual LU decomposition were used instead of the LDPU here, the complexity could be brought to \(\mathcal{O}(dnr^\omega)\) where \(\omega\) is the currently best-known power such that matrix-matrix multiplication can be done in \(\mathcal{O}(n^\omega)\). It is currently not known whether the LDPU decomposition can be done in \(\mathcal{O}(n^\omega)\) time, as is the case for the PLU decomposition. This is left as an open question. In the case of \(\mathbb{F}_2\), it might be possible to accelerate the average-case complexity by adapting the so-called “method of 4 Russians” [32] to get the LDPU decomposition. This is again left as an open question.

As a direct consequence of the last part of the theorem:

In the worst cases, the size of a tensor train decomposition will remain exponential with respect to the order of the tensor. This is akin to BDDs, where, in the worst case, the number of edges is exponential in the number of variables \(n\). In the latter case, however, the number of edges is no more than twice the number of valuations of its variables (i.e. \(2\cdot2^n\) for \(n\) variables). We may wonder whether we can similarly upper-bound the size of a reduced tensor train decomposition, with respect to the total number of entries in the tensor, i.e. \(n_1\cdot...\cdot n_d\). It turns out that here as well, there is only a small constant multiplicative overhead in the worst case:

6 Conclusion↩︎

We have presented here a normal form and a reduction strategy that minimises tensor trains over arbitrary fields in polynomial time. This strategy is based on the LDPU decomposition defined above, and for which we gave a naive implementation. This implementation is akin to the usual Gaussian elimination, and hence likely suffers from two drawbacks: i/ numerical instability when working with floating point arithmetic, and ii/ exponential bit complexity when working with exact representations of values in infinite fields. Those two drawbacks disappear when working with finite fields, such as \(\mathbb{F}_2\), when one wants to generalise decision diagrams. Addressing these two issues in infinite fields is left as an open problem.

We have also given an algorithm that directly decomposes arbitrary tensors into the unique equivalent tensor train in normal form. This algorithm and the one for reduction, the main proofs of the claims, and the definition of the normal form itself, were all made much simpler thanks to the formal graphical representation of the category \(\mathbf{Mat}_F\) and manipulation of its diagrams through the equational theory. Finally, we’ve shown that we can efficiently extract the first non-zero value and corresponding index, and that in the worst case, a tensor train uses at most \(8/3\) times the space required for a naive storage, akin to the BDD case.

Acknowledgments↩︎

The author acknowledges support from the PEPR integrated project EPiQ ANR-22-PETQ-0007 part of Plan France 2030, the ANR projects TaQC ANR-22-CE47-0012 and HQI ANR-22-PNCQ-0002.

7 Algorithms↩︎

8 Details of 3↩︎

We apply the algorithm on : \[\vec{T}^\intercal = \left(\begin{array}{rrrrrrrrrrrr} 0 & 0 & 1 & -1 & 1 & 2 & 2 & -1 & -1 & 1 & 0 & 2 \end{array}\right), \quad [2,2,3],\quad [~],\quad\begin{pmatrix}1\end{pmatrix}\] The different steps are then:

1. \(\begin{aligned}[t] &\operatorname{LDPU}\left[\begin{pmatrix} 0 & 0 & 1 & -1 & 1 & 2 \\ 2 & -1 & -1 & 1 & 0 & 2 \end{pmatrix}\right]\\ &= \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix} \overbracket{ \overbracket{\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix}}^{P'} \begin{pmatrix} 1 & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & 0 & 1 \\ 0 & 0 & 1 & -1 & 1 & 2 \end{pmatrix} }^{P'U} \end{aligned}\) \[c_1[0]=\begin{pmatrix}0&1\end{pmatrix}\qquad c_1[1]=\begin{pmatrix}2&0\end{pmatrix}\]

2. \(\begin{aligned}[t] &\operatorname{LDPU}\left[\begin{pmatrix} 0 & 0 & 1 \\ -1 & 1 & 2 \\ 1 & -\frac{1}{2} & -\frac{1}{2} \\ \frac{1}{2} & 0 & 1 \end{pmatrix}\right]\\ &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -1 & 1 \\ \frac{1}{2} & -\frac{1}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & \frac{1}{2} \end{pmatrix} \overbracket{ \overbracket{\begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}}^{P'} \begin{pmatrix} 1 & -1 & -2 \\ 0 & 1 & 3 \\ 0 & 0 & 1 \end{pmatrix} }^{P'U} \end{aligned}\) \[c_2[0]=\begin{pmatrix}1&\frac{1}{2}&0\\0&0&1\end{pmatrix}\qquad c_2[1]=\begin{pmatrix}\frac{1}{2}&\frac{1}{2}&\frac{1}{2}\\-1&0&0\end{pmatrix}\]

3. \(\begin{aligned}[t] &\operatorname{LDPU}\left[\begin{pmatrix} 0 & 0 & 1 & 1 & -1 & -2 & 0 & 1 & 3 \end{pmatrix}^\intercal\right]\\ &= \begin{pmatrix} 0 & 0 & 1 & 1 & -1 & -2 & 0 & 1 & 3 \end{pmatrix}^\intercal \begin{pmatrix} 1 \end{pmatrix} \overbracket{ \overbracket{\begin{pmatrix} 1 \end{pmatrix}}^{P'} \begin{pmatrix} 1 \end{pmatrix} }^{P'U} \end{aligned}\) \[c_3[0]=\begin{pmatrix}1\\0\\0\end{pmatrix}\qquad c_3[1]=\begin{pmatrix}-1\\1\\0\end{pmatrix}\qquad c_3[2]=\begin{pmatrix}-2\\3\\1\end{pmatrix}\]

For a quick sanity check, one can check that we can indeed compute the last \(2\) values of \(T\), indexed respectively by \((1,1,1)\) and \((1,1,2)\), as: \[c_1[1]c_2[1]c_3[1]=\begin{pmatrix}0\end{pmatrix}\qquad \text{and}\qquad c_1[1]c_2[1]c_3[2]=\begin{pmatrix}2\end{pmatrix}\]

9 Proofs↩︎

References↩︎