Definition 2.1. Let d ∈ N . Then n = ( n 1 , … , n d ) ∈ N d is called a dimension tuple of order d and x ∈ R n 1 × ⋯ × n d = : R n is called a tensor of order d and dimension n . Let N n = { 1 , … , n } then a tuple ( l 1 , … , l d ) ∈ N n 1 × ⋯ × N n d = : N n is called a multi-index and the corresponding entry of x is denoted by x (l ₁, …, l _d ). The positions 1, …, d of the indices l ₁, …, l _d in the expression x (l ₁, …, l _d ) are called modes of x.

To define further operations on tensors it is often useful to associate each mode with a symbolic index.

Definition 2.2. A symbolic index i of dimension n is a placeholder for an arbitrary but fixed natural number between 1 and n. For a dimension tuple n of order d and a tensor x ∈ R n we may write x (i ₁, …, i _d ) and tacitly assume that i _k are indices of dimension n _k for each k = 1, …, d. When standing for itself this notation means x ( i 1 , … , i d ) = x ∈ R n and may be used to slice the tensor x i 1 , l 2 , … , l d ∈ R n 1 where l k ∈ N n k are fixed indices for all k = 2, …, d. For any dimension tuple n of order d we define the symbolic multi-index i ⁿ = (i ₁, …, i _d ) of dimension n where i _k is a symbolic index of dimension n _k for all k = 1, …, d.

Remark 2.3. We use the letters i and j (with appropriate subscripts) for symbolic indices while reserving the letters k, l and m for ordinary indices.

Example 2.4. Let x be an order 2 tensor with mode dimensions n ₁ and n ₂, i.e. an n ₁-by-n ₂ matrix. Then x (ℓ ₁, j) denotes the ℓ ₁-th row of x and x (i, ℓ ₂) denotes the ℓ ₂-th column of x.

Inspired by Einstein notation we use the concept of symbolic indices to define different operations on tensors.

Definition 2.5. Let i ₁ and i ₂ be (symbolic) indices of dimension n ₁ and n ₂, respectively and let φ be a bijection φ : N n 1 × N n 2 → N n 1 n 2 .

We then define the product of indices with respect to φ as j = φ(i ₁, i ₂) where j is a (symbolic) index of dimension n ₁ n ₂. In most cases the choice of bijection is not important and we will write i ₁ ⋅ i ₂≔φ(i ₁, i ₂) for an arbitrary but fixed bijection φ. For a tensor x of dimension (n ₁, n ₂) the expression y i 1 ⋅ i 2 = x i 1 , i 2 defines the tensor y of dimension (n ₁ n ₂) while the expression x i 1 , i 2 = y i 1 ⋅ i 2 defines x ∈ R n 1 × n 2 from y ∈ R n 1 n 2 .

Definition 2.6. Consider the tensors x ∈ R n 1 × a × n 2 and y ∈ R n 3 × b × n 4 . Then the expression z i n 1 , i n 2 , j 1 , j 2 , i n 3 , i n 4 = x i n 1 , j 1 , i n 2 ⋅ y i n 3 , j 2 , i n 4 (1) defines the tensor z ∈ R n 1 × n 2 × a × b × n 3 × n 4 in the obvious way. Similary, for a = b the expression z i n 1 , i n 2 , j , i n 3 , i n 4 = x i n 1 , j , i n 2 ⋅ y i n 3 , j , i n 4 (2) defines the tensor z ∈ R n 1 × n 2 × a × n 3 × n 4 . Finally, also for a = b the expression z i n 1 , i n 2 , i n 3 , i n 4 = x i n 1 , j , i n 2 ⋅ y i n 3 , j , i n 4 (3) defines the tensor z ∈ R n 1 × n 2 × n 3 × n 4 as z i n 1 , i n 2 , i n 3 , i n 4 = ∑ k = 1 a x i n 1 , k , i n 2 ⋅ y i n 3 , k , i n 4 .

We choose this description mainly because of its simplicity and how it relates to the implementation of these operations in the numeric libraries numpy [33] and xerus [34].

This section will introduce the concept of tensor networks [35] and a graphical notation for certain operations which will simplify working with these structures. To this end we reformulate the operations introduced in the last section in terms of nodes, edges and half-edges.

Definition 2.7. For a dimension tuple n of order d and a tensor x ∈ R n the graphical representation of x is given by.where the node represents the tensor and the half-edges represent the d different modes of the tensor illustrated by the symbolic indices i ₁, …, i _d .

With this definition we can write the reshapings of Defintion 2.5 simply as

and also simplify the binary operations of Definition 2.6.

Definition 2.8. Let x ∈ R n 1 × a × n 2 and y ∈ R n 3 × b × n 4 be two tensors. Then Operation Eq. 1 is represented by

and defines z ∈ R ⋯ × a × b × ⋯ . For a = b Operation Eq. 2 is represented by

and defines z ∈ R ⋯ × a × ⋯ and Operation Eq. 3 defines z ∈ R ⋯ × ⋯ by.

With these definitions we can compose entire networks of multiple tensors which are called tensor networks.

A prominent example of a tensor network is the tensor train (TT) [19,36], which is the main tensor network used throughout this work. This network is discussed in the following subsection.

Definition 2.9. Let n be an dimensional tuple of order-d. The TT format decomposes an order d tensor x ∈ R n into d component tensors x k ∈ R r k − 1 × n k × r k for k = 1, …, d with r ₀ = r _d = 1. This can be written in tensor network formula notation as x i 1 , ⋯ , i d = x 1 i 1 , j 1 ⋅ x 2 j 1 , i 2 , j 2 ⋯ x d j d − 1 , i d .

The tuple (r ₁, …, r _d−1) is called the representation rank of this representation.

In graphical notation it looks like this.

Remark 2.10. Note that this representation is not unique. For any pair of matrices (A, B) that satisfies AB = Id we can replace x _k by x _k (i ₁, i ₂, j) ⋅ A (j, i ₃) and x _k+1 by B (i ₁, j) ⋅ x (j, i ₂, i ₃) without changing the tensor x.

The representation rank of x is therefore dependent on the specific representation of x as a TT, hence the name. Analogous to the concept of matrix rank we can define a minimal necessary rank that is required to represent a tensor x in the TT format.

Definition 2.11. The tensor train rank of a tensor x ∈ R n with tensor train components x 1 ∈ R n 1 × r 1 , x k ∈ R r k − 1 × n k × r k for k = 2, …, d − 1 and x d ∈ R r d − 1 × n d is the set TT - rank x = r 1 , … , r d of minimal r _k ’s such that the x _k compose x.

In [[22], Theorem 1a] it is shown that the TT-rank can be computed by simple matrix operations. Namely, r _k can be computed by joining the first k indices and the remaining d − k indices and computing the rank of the resulting matrix. At last, we need to introduce the concept of left and right orthogonality for the tensor train format.

Definition 2.12. Let x ∈ R m × n be a tensor of order d + 1. We call x left orthogonal if x i m , j 1 ⋅ x i m , j 2 = Id j 1 , j 2 .

Similarly, we call a tensor x ∈ R m × n of order d + 1 right orthogonal if x i 1 , j n ⋅ x i 2 , j n = Id i 1 , i 2 .

A tensor train is left orthogonal if all component tensors x ₁, …, x _d−1 are left orthogonal. It is right orthogonal if all component tensors x ₂, …, x _d are right orthogonal.

Lemma 2.1 [36]. For every tensor x ∈ R n of order d we can find left and right orthogonal decompositions.

For technical purposes it is also useful to define the so-called interface tensors, which are based on left and right orthogonal decompositions.

Definition 2.13. Let x be a tensor train of order d with rank tuple r .

For every k = 1, …, d and ℓ = 1, …, r _k , the ℓ-th left interface vector is given by τ k , ℓ ≤ ( x ) i 1 , i 2 , … , i k = x 1 i 1 , j 1 … x k j k − 1 , i k , ℓ where x is assumed to be left orthogonal. The ℓ-th right interface vector is given by τ k + 1 , ℓ ≥ ( x ) i k + 1 , … , i d = x k + 1 ℓ , i k + 1 , j k + 1 … x d j d − 1 , i d where x is assumed to be right orthogonal.

In this section we specify the setup for our method and define the majority of the different sets of polynomials that are used. We start by defining dictionaries of one dimensional functions which we then use to construct the different sets of high-dimensional functions.

Definition 2.14. Let p ∈ N be given. A function dictionary of size p is a vector valued function Ψ = ( Ψ 1 , … , Ψ p ) : R → R p .

Example 2.15. Two simple examples of a function dictionary that we use in this work are given by the monomial basis of dimension p, i.e. Ψ monomial x = 1 x x 2 … x p − 1 T (4) and by the basis of the first p Legendre polynomials, i.e. Ψ Legendre x = 1 x 1 2 3 x 2 − 1 1 2 5 x 3 − 3 x … T . (5)

Using function dictionaries we can define the following high-dimensional space of multivariate functions. Let Ψ be a function dictionary of size p ∈ N . The d-th order product space that corresponds to the function dictionary Ψ is the linear span V p d ≔ ⊗ k = 1 d Ψ m k : m ∈ N p d . (6)

This means that every function u ∈ V p d can be written as u x 1 , … , x d = c i 1 , … , i d ∏ k = 1 d Ψ x k i k (7) with a coefficient tensor c ∈ R p where p = (p, …, p) is a dimension tuple of order d. Note that equation Eq. 7 uses the index notation from Definition 2.6 with arbitrary but fixed x _k ’s. Since R p is an intractably large space, it makes sense for numerical purposes to consider the subset T r V p d ≔ u ∈ V p d : TT - rank c ≤ r (8) where the TT rank of the coefficient is bounded. Every u ∈ T r ( V p d ) can thus be represented graphically as (9) where the C _k ’s are the components of the tensor train representation of the coefficient tensor c ∈ R p of u ∈ V p d .

Remark 2.16. In this way every tensor c ∈ R p (in the tensor train format) corresponds one to one to a function u ∈ V p d .

An important subspace of V p d is the space of homogeneous polynomials. For the purpose of this paper we define the subspace of homogeneous polynomials of degree g as the space W g d ≔ ⊗ k = 1 d Ψ m k : m ∈ N p d and ∑ k = 1 d m k = d + g , (10) where again ⟨•⟩ is the linear span. From this definition it is easy to see that a homogeneous polynomial of degree g can be represented as an element of V p d where the coefficient tensor c satisfies c m 1 , … , m d = 0 , i f ∑ k = 1 d m k ≠ d + g .

In Theoretical Foundation we will introduce an efficient representation of such coefficient tensors c in a block sparse tensor format.

Using W g d we can also define the space of polynomials of degree at most g by S g d = ⊕ g ̃ = 0 g W g ̃ d . (11)

Based on this characterization we will define a block-sparse tensor train version of this space in Theoretical Foundation.

In algebraic geometry the space W g d is considered classically only for the dictionary Ψ_monomial of monomials and is typically parameterized by a symmetric tensor u x = B i 1 , ··· , i g ⋅ x i 1 ··· x i g , x ∈ R d (12) where d = (d, …, d) is a dimension tuple of order g and B ∈ R d satisfies B (m ₁, …, m _g ) = B (σ(m ₁, …, m _g )) for every permutation σ in the symmetric group S _g . We conclude this section by showing how the representation Eq. 7 can be calculated from the symmetric tensor representation Eq. 12, and vice versa. By equating coefficients we find that for every ( m 1 , … , m d ) ∈ N p d either m ₁ + ⋯ + m _d ≠ d + g and c (m ₁, …, m _d ) = 0 or c m 1 , … , m d = ∑ σ n : σ ∈ S g B σ n 1 , … , n g where n 1 , … , n g = ( 1 , … , 1 ︸ m 1 − 1  times , 2 , … , 2 ︸ m 2 − 1  times , … , ) ∈ N d g .

Since B is symmetric the sum simplifies to ∑ σ n : σ ∈ S g B σ n 1 , … , n g = g m 1 − 1 , … , m d − 1 B n 1 , … , n g .

From this follows that for ( n 1 , … , n g ) ∈ N d g B n 1 , … , n g = 1 g m 1 − 1 , … , m d − 1 c m 1 , … , m d where m k = 1 + ∑ ℓ = 1 g δ k , n ℓ for all k = 1 , … , d and δ_k,ℓ denotes the Kronecker delta. This demonstrates how our approach can alleviate the difficulties that arise when symmetric tensors are represented in the hierarchical tucker format [37] in a very simple fashion.

Let in the following V p d be the product space of a function dictionary Ψ such that V p d ⊆ L 2 ( Ω ) . Consider a high-dimensional function f ∈ L ₂(Ω) on some domain Ω ⊂ R d and assume that the point-wise evaluation f(x) is well-defined for x ∈ Ω. In practice it is often possible to choose Ω as a product domain Ω = Ω₁ × Ω₂ ×⋯ Ω _d by extending f accordingly. To find the best approximation u _W of f in the space W ⊆ V p d we then need to solve the problem u W = a r g m i n u ∈ W ‖ f − u ‖ L 2 Ω 2 . (13)

A practical problem that often arises when computing u _W is that computing the L ₂(Ω)-norm is intractable for large d. Instead of using classical quadrature rules one often resorts to a Monte Carlo estimation of the high-dimensional integral. This means one draws M random samples { x ( m ) } m = 1 , … , M from Ω and estimates ‖ f − u ‖ L 2 Ω 2 ≈ 1 M ∑ m = 1 M ‖ f x m − u x m ‖ F 2 , where ‖⋅‖ _F is the Frobenius norm. With this approximation we can define an empirical version of u _W as u W , M = a r g m i n u ∈ W 1 M ∑ m = 1 M ‖ f x m − u x m ‖ F 2 . (14)

For a linear space W, computing u _W,M amounts to solving a linear system and does not pose an algorithmic problem. We use the remainder of this section to comment on the minimization problem Eq. 14 when a set of tensor trains is used instead.

Given samples ( x ( m ) ) m = 1 , … , M we can evaluate u ∈ V p d for each x ( m ) = ( x 1 ( m ) , … , x d ( m ) ) using Eq. 7 If the coefficient tensor c of u can be represented in the TT format then we can use Eq. 9 to perform this evaluation efficiently for all samples ( x ( m ) ) m = 1 , … , M at once. For this we introduce for each k = 1, …, d the matrix Ξ k = Ψ x k 1 … Ψ x k M ∈ R p × M . (15)

Then the M-dimensional vector of evaluations of u at all given sample points is given by.

where we use Operation Eq. 2 to join the different M-dimensional indices.

The alternating least-squares algorithm cyclically updates each component tensor C _k by minimizing the residual corresponding to this contraction. To formalize this we define the operator Φ k ∈ R M × r k − 1 × n k × r k as . (16)

Then the update for C _k is given by a minimal residual solution of the linear system C k = a r g m i n C ∈ R r k − 1 × n k × r k ‖ Φ k j , i 1 , i 2 , i 3 ⋅ C i 1 , i 2 , i 3 − F j ‖ F 2 where F(m)≔y ^(m)≔f (x ^(m)) and i ₁, i ₂, i ₃, j are symbolic indices of dimensions r _k−1, n _k , r _k , M, respectively. The particular algorithm that is used for this minimization may be adapted to the problem at hand. These contractions are the basis for our algorithms in Method Description. We refer to [19] for more details on the ALS algorithm.

Note that it is possible to reuse parts of the contractions in Φ _k through so called stacks. In this way not the entire contraction has to be computed for every k. The dashed boxes mark the parts of the contraction that can be reused. Details on that can be found in [38].

The accuracy of the solution u _W,M of Eq. 14 in relation to u _W is subject to tremendous interest on the part of the mathematics community. Two particular papers that consider this problem are [32,39]. While the former provides sharper error bounds for the case of linear ansatz spaces the latter generalizes the work and is applicable to tensor network spaces. We now recall the relevant result for convenience.

Proposition 3.1. For any set W ⊆ L ²(Ω) ∩ L ^∞ (Ω), define the variation constant K W ≔ sup v ∈ W \ 0 ‖ v ‖ L ∞ Ω 2 ‖ v ‖ L 2 Ω 2 .

Let δ ∈ (0, 2^−1/2). If W is a subset of a finite dimensional linear space and k≔ max{K ({f − u _W }), K ({u _W }− W)} < ∞ it holds that P ‖ f − u W , M ‖ L 2 Ω ≤ 3 + 4 δ ‖ f − u W ‖ L 2 Ω ≥ 1 − q where q decreases exponentially with a rate of ln ( q ) ∈ O ( − M δ 2 k − 2 ) .

Proof. Since k < ∞, Theorems 2.7 and 2.12 in [32] ensure that ‖ f − u W , M ‖ L 2 Ω ≤ 1 + 2 1 + δ 1 − δ ‖ f − u W ‖ L 2 Ω . holds with a probability of at least 1 − 2 C ⁡ exp ( − 1 2 M δ 2 k − 2 ) . The constant C is independent of M and, since W is a subset of a finite dimensional linear space, depends only polynomially on δ and k ⁻¹. For δ ∈ (0, 2^−1/2) it holds that 1 + δ 1 − δ ≤ 1 + 2 δ . This concludes the proof.

Note that the value of k depends only on f and on the set W but not on the particular choice of representation of W. However, the variation constant of spaces like V p d still depends on the underlying dictionary Ψ. Although the proposition indicates that a low value of k is necessary to achieve a fast convergence, the tensor product spaces V p d considered thus far does not exhibit a small variation constant. The consequence of Proposition 3.1 is that elements of this space are hard to learn in general and may require an infeasible number of samples. To see this consider Ω = [−1,1] ^d and the function dictionary Ψ_Legendre of Legendre polynomials Eq. 5. Let L ⊆ N p d and define P ℓ ( x ) ≔ ∏ k = 1 d 2 ℓ k − 1 ( Ψ Legendre ( x k ) ) ℓ k for all ℓ ∈ L. Then, { P ℓ } ℓ ∈ L is an L ²-orthonormal basis for the linear subspace V ≔ ⟨ P ℓ : ℓ ∈ L ⟩ ⊆ V p d and one can show that K V = sup x ∈ Ω ∑ ℓ ∈ L P ℓ x 2 = ∑ ℓ ∈ L ∏ k = 1 d 2 ℓ k − 1 , (17) by using techniques from [[32], Sample Complexity for Polynomials] and the fact that each P _ℓ attains its maximum at 1. If L = N p d , we can interchange the sum and product in Eq. 17 and can conclude that K ( V p d ) = p 2 d . This means that we have to restrict the space V p d to obtain an admissible variation constant. We propose to use the space W g d of homogeneous polynomials of degree g. Employing Eq. 17 with L = { ℓ : | ℓ | = d + g} we obtain the upper bound K W g d ≤ d − 1 + g d − 1 max | ℓ | = d + g ∏ k = 1 d 2 ℓ k − 1 ≤ d − 1 + g d − 1 2 ⌊ g d ⌋ + 3 g mod d 2 ⌊ g d ⌋ + 1 d − g mod d where the maximum is estimated by observing that (2 (ℓ ₁ + 1) − 1) (2ℓ ₂ − 1) ≤ (2ℓ ₁ − 1) (2 (ℓ ₂ + 1) − 1) ⇔ ℓ ₂ ≤ ℓ ₁. For g ≤ d this results in the simplified bound K ( W g d ) ≤ ( 3 e d − 1 + g g ) g , where e is the Euler number. This improves the variation constant substantially compared to the bound K ( V p d ) ≤ p 2 d , when g ≪ d. A similar bound for the dictionary of monomials Ψ_monomial is more involved but can theoretically be computed in the same way.

In this work, we focus on the case where the samples are drawn according to a probability measure on Ω. This however is not a necessity and it is indeed beneficial to draw the samples from an adapted sampling measure. Doing so, the theory in [32] ensures that K(V) = dim(V) for all linear spaces V — independent of the underlying dictionary Ψ. This in turn leads to the bounds K ( V p d ) = p d and K ( W g d ) = ( d − 1 + g d − 1 ) ≤ ( e d − 1 + g g ) g for g ≤ d. These optimally weighted least-squares methods however, are not the focus of this work and we refer the interested reader to the works [39,40].

Now that we have seen that it is advantageous to restrict ourselves to the space W g d we need to find a way to do so without loosing the advantages of the tensor train format. In [29] it was rediscovered from the physics community (see [28]) that if a tensor train is an eigenvector of certain Laplace-like operators it admits a block sparse structure. This means for a tensor train c the components C _k have zero blocks. Furthermore, this block sparse structure is preserved under key operations, like e.g. the TT-SVD. One possible operator which introduces such a structure is the Laplace-like operator L = ∑ k = 1 d ⊗ ℓ = 1 k − 1 I p ⊗ diag 0,1 , … , p − 1 ⊗ ⊗ ℓ = k + 1 d I p . (18)

This is the operator mentioned in the introduction encoding a quantum symmetry. In the context of quantum mechanics this operator is known as the bosonic particle number operator but we simply call it the degree operator. The reason for this is that for the function dictionary of monomials Ψ_monomial the eigenspaces of L for eigenvalue g are associated with homogeneous polynomials of degreee g. Simply put, if the coefficient tensor c for the multivariate polynomial u ∈ V p d is an eigenvector of L with eigenvalue g, then u is homogeneous and the degree of u is g. In general there are polynomials in V p d with degree up to (p − 1)d. To state the results on the block-sparse representation of the coefficient tensor we need the partial operators L k ≤ = ∑ m = 1 k ⊗ ℓ = 1 m − 1 I p ⊗ diag 0,1 , … , p − 1 ⊗ ⊗ ℓ = m + 1 k I p L k + 1 ≥ = ∑ m = k + 1 d ⊗ ℓ = k + 1 m − 1 I p ⊗ diag 0,1 , … , p − 1 ⊗ ⊗ ℓ = m + 1 d I p , for which we have L = L k ≤ ⊗ ⊗ ℓ = k + 1 d I p + ⊗ l = 1 k I p ⊗ L k + 1 ≥ .

In the following we adopt the notation x = Lc to abbreviate the equation x i 1 , … , i d = L i 1 , … , i d , j 1 , … , j d c j 1 , … , j d where L is a tensor operator acting on a tensor c with result x.

Recall that by Remark 2.16 every TT corresponds to a polynomial by multiplying function dictionaries onto the cores. This means that for every ℓ = 1, …, r the TT τ k , ℓ ≤ ( c ) corresponds to a polynomial in the variables x ₁, …, x _k and the TT τ k + 1 , ℓ ≥ ( c ) corresponds to a polynomial in the variables x _k+1, …, x _d . In general these polynomials are not homogeneous, i.e. they are not eigenvectors of the degree operators L k ≤ and L k + 1 ≥ . But since TTs are not uniquely defined (cf. Remark 2.10) it is possible to find transformations of the component tensors C _k and C _k+1 that do not change the tensor c or the rank r but result in a representation where each τ k , ℓ ≤ ( c ) and each τ k + 1 , ℓ ≥ ( c ) correspond to a homogeneous polynomial. Thus, if c represents a homogeneous polynomial of degree g and τ k , ℓ ≤ ( c ) is homogeneous with deg ( τ k , ℓ ≤ ( c ) ) = g ̃ then τ k + 1 , ℓ ≥ ( c ) must be homogeneous with deg ( τ k , ℓ ≥ ( c ) ) = g − g ̃ .

This is put rigorously in the first assertion in the subsequent Theorem 3.2. There S k , g ̃ contains all the indices ℓ for which the reduced basis polynomials satisfy deg ( τ k , ℓ ≤ ( c ) ) = g ̃ . Equivalently, it groups the basis functions τ k + 1 , ℓ ≥ ( c ) into functions of order g − g ̃ . The second assertion in Theorem 3.2 states that we can only obtain a homogeneous polynomial of degree g ̃ + m in the variables x ₁, …, x _k by multiplying a homogeneous polynomial of degree g ̃ in the variables x ₁, …, x _k−1 with a univariate polynomial of degree m in the variable x _k . This provides a constructive argument for the proof and can be used to ensure block-sparsity in the implementation. Note that this condition forces entire blocks in the component tensor C _k in equation (20) to be zero and thus decreases the degrees of freedom.

Theorem 3.2 [[29], Theorem 3.2]. Let p = (p, …, p) be a dimension tuple of size d and c ∈ R p \ { 0 } , be a tensor train of rank r = (r ₁, …, r _d−1). Then Lc = gc if and only if c has a representation with component tensors C k ∈ R r k − 1 × p × r k that satisfies the following two properties.

1. For all g ̃ ∈ { 0,1 , … , g } there exist S k , g ̃ ⊆ { 1 , … , r k } such that the left and right unfoldings satsify L k ≤ τ k , ℓ ≤ ( c ) = g ̃ τ k , ℓ ≤ ( c ) L k + 1 ≥ τ k + 1 , ℓ ≥ ( c ) = g − g ̃ τ k + 1 , ℓ ≥ ( c ) (19) for ℓ ∈ S k , g ̃ .

2. The component tensors satisfy a block structure in the sets S k , g ̃ for m = 1, …, p C k ℓ 1 , m , ℓ 2 ≠ 0 ⇒ ∃ 0 ≤ g ̃ ≤ g − m − 1 : ℓ 1 ∈ S k − 1 , g ̃ ∧ ℓ 2 ∈ S k , g ̃ + m − 1 (20) where we set S 0,0 = S d , g = { 1 } .

Note that this generalizes to other dictionaries and is not restricted to monomials.

Although, block sparsity also appears for g + 1 ≠ p we restrict ourselves to the case g + 1 = p in this work. Note that then the eigenspace of L for the eigenvalue g has a dimension equal to the dimension of the space of homogeneous polynomials, namely ( d − 1 + g d − 1 ) . Defining ρ k , g ̃ ≔ | S k , g ̃ | , we can derive the following rank bounds.

Lemma 3.3 [[29], Lemma 3.6]. Let p = (p, …, p) be a dimension tuple of size d and c ∈ R p \ { 0 } , with Lc = gc . Assume that g + 1 = p then the block sizes ρ k , g ̃ from Theorem 3.2 are bounded by ρ k , g ̃ ≤ min k + g ̃ − 1 k − 1 , d − k + g − g ̃ − 1 d − 1 (21) for all k = 1, …, d − 1 and g ̃ = 0 , … , g and ρ _k,0 = ρ _k,g = 1.

The proof of this lemma is based on a simple combinatorial argument. For every k consider the size of the groups ρ k − 1 , g ̄ for g ̄ ≤ g ̃ . Then ρ k , g ̄ can not exceed the sum of these sizes. Similarly, ρ k , g ̄ can not exceed ∑ g ̄ ≤ g ̃ ρ k + 1 , g ̄ . Solving these recurrence relations yields the bound.

Example 3.1 (Block Sparsity). Let p = 4 and g = 3 be given and let c be a tensor train such that Lc = gc. Then for k = 2, …, d − 1 the component tensors C _k of c exhibit the following block sparsity (up to permutation). For indices i of order r _k−1 and j of order r _k C k i , 1 , j = * 0 0 0 0 * 0 0 0 0 * 0 0 0 0 * C k i , 2 , j = 0 * 0 0 0 0 * 0 0 0 0 * 0 0 0 0 C k i , 3 , j = 0 0 * 0 0 0 0 * 0 0 0 0 0 0 0 0 C k i , 4 , j = 0 0 0 * 0 0 0 0 0 0 0 0 0 0 0 0 .

This block structure results from sorting the indices i and j in such a way that max S k , g ̃ + 1 = min S k , g ̃ + 1 for every g ̃ .

The maximal block sizes ρ k , g ̃ for k = 1, …, d − 1 are given by ρ k , 0 = 1 , ρ k , 1 = min k , d − k , ρ k , 2 = min k , d − k , ρ k , 3 = 1 .

As one can see by Lemma 3.3 the block sizes ρ k , g ̃ can still be quite high.

The expressive power of tensor train parametrizations can be understood by different concepts, such as locality or self similarity. We use the remainder of this section to provide d-independent rank bounds in the context of locality.

Definition 3.2. Let u ∈ W g d be a homogeneous polynomial and B be the symmetric coefficient tensor introduced in Parametrizing homogeneous polynomials by symmetric tensors We say that u has a variable locality of K _loc if B (ℓ ₁, …, ℓ _g ) = 0 for all ( ℓ 1 , … , ℓ g ) ∈ N d g with max | ℓ m 1 − ℓ m 2 | : m 1 , m 2 = 1 , … , g > K loc .

Example 3.3. Let u be a homogeneous polynomial of degree 2 with variable locality K _loc. Then the symmetric matrix B (cf. Eq. 12) is K _loc-banded. For K _loc = 0 this means that B is diagonal and that u takes the form u x = ∑ ℓ = 1 d B ℓ ℓ x ℓ 2 .

This shows that variable locality removes mixed terms.

Remark 3.4. The locality condition in the following Theorem 3.4 is a sufficient, but in no way necessary, condition for a low rank. But since locality is a prominent feature of many physical phenomena, this condition allows us to identify an entire class of highly relevant functions which can be approximated very efficiently.

Consider, for example, a many-body system in one dimension, where each body is described by position and velocity coordinates. If the influence of neighboring bodies is much higher than the influence of more distant ones, the coefficients of the polynomial parts that depend on multiple variables often can be neglected. The forces in this system then exhibit a locality structure. An example of this is given in equation Eq. 6 in [3], where this structure is exhibited by the force that acts on the bodies. A similar structure also appears in the microscopic traffic models in Notation of [41].

Another example is given by the polynomial chaos expansion of the stochastic process X t ξ 1 , … , ξ t ≔ c i 1 , … , i t ∏ k = 1 t Ψ ξ k i k for t ∈ N , where Ψ is the function dictionary of Hermite polynomials. In many applications, it is justified to assume that the magnitude of the covariance Cov ( X t 1 , X t 2 ) decays with the distance of the indices |t ₁ − t ₂| If the covariance decays fast enough, the coefficient tensor exhibits approximate locality, i.e. it can be well approximated by a coefficient tensor that satisfies the locality condition. Examples of this are Gaussian processes with a Matérn kernel [42,43] or Markov processes.

Theorem 3.4. Let p = (p, …, p) be a dimension tuple of size d and c ∈ R p \ { 0 } correspond to a homogeneous polynomial of degree g + 1 = p (i.e. Lc = gc ) with variable locality K _loc . Then the block sizes ρ k , g ̃ are bounded by ρ k , g ̃ ≤ ∑ ℓ = 1 K loc min { K loc − ℓ + 1 + g ̃ − 2 K loc − ℓ , ℓ + g − g ̃ − 2 ℓ − 1 } (22) for all k = 1, …, d − 1 and g ̃ = 1 , … , g − 1 as well as ρ _k,0 = ρ _k,g = 1.

Proof. For fixed g > 0 and a fixed component C _k recall that for each l the tensor τ k , l ≤ ( c ) corresponds to a reduced basis function v _l in the variables x ₁, …, x _k and that for each l the tensor τ k + 1 , l ≥ ( c ) corresponds to a reduced basis function w _l in the variables x _k+1, …, x _d . Further recall that the sets S k , g ̃ group these v _l and w _l . For all l ∈ S k , g ̃ it holds that deg ( v l ) = g ̃ and deg ( w l ) = g − g ̃ . For g ̃ = 0 and g ̃ = g we know from Lemma 3.3 that ρ k , g ̃ = 1 . Now fix any 0 < g ̃ < g and arrange all the polynomials v _l of degree g ̃ in a vector v and all polynomials w _l of degree g − g ̃ in a vector w. Then every polynomial of the form v ^⊺ Qw for some matrix Q satisfies the degree constraint and the maximal possible rank of Q provides an upper bound for the block size ρ k , g ̃ . However, due to the locality constraint we know that certain entries of Q have to be zero. We denote a variable of a polynomial as inactive if the polynomial is constant with respect to changes in this variable and active otherwise. Assume that the polynomials in v are ordered (ascendingly) according to the smallest index of their active variables and that the polynomials in w are ordered (ascendingly) according to the largest index of their active variables. With this ordering Q takes the form Q = 0 Q 1 * Q 2 * * Q 3 ⋮ ⋮ ⋮ ⋱ * * * ⋯ Q K l o c .

This means that for ℓ = 1, …, K _loc each block Q _ℓ matches a polynomial v _l of degree g ̃ in the variables x k − K loc + ℓ , … , x k with a polynomial w _l of degree g − g ̃ in the variables x _k+1, …, x _k+ℓ .

Observe that the number of rows in Q _ℓ decreases while the number columns increases with ℓ. This means that we can subdivide Q as Q = 0 0 0 Q C 0 0 * Q R 0 , where Q _C contains the blocks Q _ℓ with more rows than columns (i.e. full column rank) and Q _R contains the blocks Q _ℓ with more columns than rows (i.e. full row rank). So Q _C is a tall-and-skinny matrix while Q _R is a short-and-wide matrix and the rank for general Q is bounded by the sum over the column sizes of the Q _ℓ in Q _C plus the sum over the row sizes of the Q _ℓ in Q _R i.e. rank Q = ∑ ℓ = 1 K loc rank Q ℓ .

To conclude the proof it remains to compute the row and column sizes of Q _ℓ . Recall that the number of rows of Q _ℓ equals the number of polynomials u of degree g ̃ in the variables x k − K loc + ℓ , … , x k that can be represented as u ( x k − K loc + ℓ , … , x k ) = x k − K loc + ℓ u ̃ ( x k − K loc + ℓ , … , x k ) . This corresponds to all possible u ̃ of degree g ̃ − 1 in the K _loc − ℓ + 1 variables x k − K loc + ℓ , … , x k . This means that #rows Q ℓ ≤ K loc − ℓ + 1 + g ̃ − 2 K loc − ℓ and a similar argument yields #columns Q ℓ ≤ ℓ + g − g ̃ − 2 ℓ − 1 .

This concludes the proof.

This lemma demonstrates how the combination of the model space W g d with a tensor network space can reduce the space complexity by incorporating locality.

Remark 3.5. The rank bound in Theorem 3.4 is only sharp for the highest possible rank. The precise bounds can be much smaller, especially for the first and last ranks, but are quite technical to write down. For this reason, we do not provide them.

One sees that the bound only depends on g and K _loc and is therefore d-independent.

Remark 3.6. The rank bounds presented in this section do not only hold for the monomial dictionary Ψ_monomial but for all polynomial dictionaries Ψ that satisfy deg(Ψ _k ) = k − 1 for all k = 1, …, p. When we speak of polynomials of degree g, we mean the space W g d = { v ∈ V p d : deg ( v ) = g } . For the dictionary of monomials Ψ_monomial the space W g d contains only homogeneous polynomials in the classical sense. However, when the basis of Legendre polynomials Ψ_Legendre is used one obtains a space in which the functions are not homogeneous in the classical sense. Note that we use polynomials since they have been applied successfully in practice, but other function dictionaries can be used as well. Also note that the theory is much more general as shown in [29] and is not restricted to the degree counting operator.

With Theorem 3.4 one sees that tensor trains are well suited to parametrize homogeneous polynomials of fixed degree where the symmetric coefficient tensor B (cf. Eq. 12) is approximately banded (see also Example 3.3). This means, that there exist an K _loc such that the error for a best approximation of B by a tensor B ̃ with variable locality K _loc is small. However, K _loc is not known precisely in practice but can only be assumed by physical understanding of the problem at hand. Therefore, we still rely on rank adaptive schemes to find appropriate rank and block sizes. Moreover, the locality property heavily depends on the ordering of the modes. This ordering can be optimized, for example, by using entropy measures for the correlation of different modes, as it is done in quantum chemistry (cf. [[44], Remark 4.2]) or by model selection methods (cf. [25,27,45,46]).