Approximation in the extended functional tensor train format

This work proposes the extended functional tensor train (EFTT) format for compressing and working with multivariate functions on tensor product domains. Our compression algorithm combines tensorized Chebyshev interpolation with a low-rank approximation algorithm that is entirely based on function evaluations. Compared to existing methods based on the functional tensor train format, our approach often reduces the required storage, sometimes considerably, while achieving the same accuracy. In particular, we reduce the number of function evaluations required to achieve a prescribed accuracy by up to over 96% compared to the algorithm from [Gorodetsky, Karaman and Marzouk, Comput. Methods Appl. Mech. Eng. , 347 (2019)] .


Introduction
Multivariate functions become notoriously difficult to approximate and work with in high dimensions.The storage and computational effort of standard approximation techniques increases exponentially with the dimension.Several techniques have been developed to potentially mitigate this so-called curse of dimensionality, including sparse grids [15] and low-rank tensor approximation techniques [46].On a functional level, low-rank tensor formats correspond to separation of variables.In particular, if f : [−1, 1] d → R can be written as a sum of separable functions, f (x 1 , . . ., x d ) = R α=1 g (1)  α (x 1 )g (2)  α (x 2 ) then the task of approximating f is replaced by approximating the R•d univariate functions g (ℓ) α : [−1, 1] → R; see, e.g., [12].It is often beneficial, especially when R is very large, to impose additional structure on (1), leading to the so-called functional tensor train (FTT) [13,38] and Tucker [28,49] formats.In practice, functions are rarely given directly in such a functional low-rank format, but they can sometimes be well approximated by it.This is exploited in quite a few applications, including the solution of time-dependent partial differential equations (PDEs) [22], uncertainty quantification [56], sensitivity analysis [9], optimal control [37] and quantum dynamic simulations [81].
The FTT format [13,38] takes the form 1,α1 (x 1 )g (2)  α1,α2 (x 2 ) for univariate functions g α ℓ−1 ,α ℓ by, e.g., a truncated series expansion, which incurs additional truncation error [38].It is important to emphasize that not every function admits a good approximation in the FTT format (2) with moderate TT ranks, especially for larger d.Generally speaking, smoothness [47] and a restricted (e.g., nearest neighbor) interaction between variables are helpful.Upper bounds for the TT ranks needed to attain a certain accuracy are derived in [13,42] for functions in Sobolev spaces.The required TT ranks can change significantly when variables are transformed [91] or the order of variables is permuted [23].For functions in periodic and mixed Sobolev spaces approximation rates can be found in [79] and [43], respectively.Much better rates can be obtained for functions with special structures such as compositional functions [7] and quantities associated with certain parametric PDEs [6].For a more abstract analysis of functional low-rank approximation we refer to [1,2,3].
To compute functional low-rank approximations in practice, one often starts from multivariate polynomial interpolation [60,92].When using a tensorized grid of interpolation nodes, the number of function evaluations and the memory needed to store the coefficients of the interpolant grow exponentially with d.This growth can be alleviated by utilizing a low-rank approximation of the coefficient tensor, which corresponds to a functional low-rank approximation [28,46].Such an approach is used for bivariate [88] and trivariate [28,49] functions by Chebfun [29] and for general multivariate functions in [13].The impact of a TT approximation of the coefficient tenor on the overall interpolation error has been studied for tensorized Lagrange interpolation with respect to the L 2 -norm in [13].In [28], the impact of arbitrary coefficient tensor approximations is analyzed with respect to the uniform norm in the context of tensorized Chebyshev interpolation of trivariate functions.A slightly different approach to obtain FTT approximations is proposed in [38], which applies a continuous version of the so-called TT-cross algorithm [67] directly to the function.All univariate functions occurring in this process are discretized using parameterizations such as basis expansions.In principle, the format allows every univariate function to be parameterized differently.In the special case that univariate interpolation with the same basis functions is used to discretize all univariate functions in the same mode, this approach is equivalent to a TT approximation of the coefficient tensor for tensorized interpolation.
The FTT format (2) requires the storage of R ℓ−1 • R ℓ univariate functions for each ℓ = 1, . . ., d.These functions (or their parametrizations) will often be linearly dependent.In this case, we can compress the format further by storing these functions in terms of linear combinations of r ℓ ≪ R ℓ−1 • R ℓ basis functions.The resulting approximation can be seen as a variant of the hierarchical Tucker format proposed in [48].When the approximation is constructed from a low-rank approximation of the coefficient tensor, we can obtain such a compressed format by following the ideas of Khoromskij [52], who suggested to use a two-level Tucker format, which is obtained by first computing a Tucker approximation [94] before approximating the core tensor in R r1×•••×r d further.Computing a TT approximation of the core tensor leads to a so-called extended TT format [27,30,79], from which we can derive a functional low-rank approximation in the extended functional tensor train (EFTT) format.This EFTT approximation corresponds to a compressed FTT approximation.
In this work, we propose a novel algorithm to efficiently compute functional low-rank approximations using the EFTT format.Our algorithm is based on first obtaining suitable factor matrices for a Tucker approximation of the coefficient tensor.The columns of these matrices are determined by applying a combination of adaptive cross approximation [10] with randomized pivoting to matricizations of the coefficient tensor.We then apply oblique projections based on discrete empirical interpolation [16] as in [28,82] to implicitly construct a suitable core tensor for the Tucker approximation.In a second step, we compute a TT approximation of this core tensor using the greedy restricted cross interpolation algorithm [78], which is a rank adaptive variant of the TT-cross algorithm [67] requiring asymptotically fever evaluations than the TT-DMRG-cross algorithm [77].The combination of the Tucker and TT approximation yields the desired functional low-rank approximation.The main advantage of our approach compared to a direct TT approximation of the coefficient tensor is that the approximated core tensor is potentially much smaller than the coefficient tensor.This reduces the storage complexity of the approximation and it significantly reduces the number of function evaluations needed in the computation of the TT approximation, which we demonstrate in our numerical experiments.Other approaches that can potentially achieve this goal have been discussed in [41,44,62,93].
Our novel algorithm is potentially useful in all applications that require to store and work with a multivariate functions.This includes uncertainty quantification, computational chemistry and physics; we refer to the surveys [40,51,53] for concrete examples.
The TT format considered in this work corresponds to a line graph tensor network topology.Other tensor network topologies have been considered in the literature for the purpose of function approximation.In [8], a variant of the TT-cross algorithm for tree tensor networks, corresponding to the hierarchical Tucker format, has been developed.In [45], strategies for choosing a good tree tensor network have been proposed.It would be interesting to compare and possible combine the developments from [45] with this work.
The remainder of this paper is structured as follows.In Section 2, we recall how functional low-rank approximations can be obtained by combining tensorized Chebyshev interpolation and low-rank approximations of the coefficient tensor.Our novel algorithm to compute approximations in the EFTT format is presented in Section 3. In our numerical experiments in Section 4, we demonstrate the advantages of our novel algorithm and compare it to a direct TT approximation of the coefficient tensor as in [13] and to the FTT algorithm in [38].
2 Functional low-rank approximation via tensorized interpolation

Tensorized Chebyshev interpolation
Given a function f : [−1, 1] d → R, let us consider a polynomial approximation of degree (n 1 , . . ., n d ) taking the form where ) is the coefficient tensor and T k (x) = cos(k cos −1 (x)) denotes the k-th Chebyshev polynomial.To be consistent with standard notation [65], we index entries of the coefficient tensor starting from 0. Matrices and tensors not related to Chebyshev approximation are indexed in the usual way, starting from 1.
For analytic functions f it is shown in [76,Lemma 7.3.3]that the error of the approximation (3) in the uniform norm decays exponentially with respect to min{n 1 , . . ., n d }.

Functional tensor train (FTT) approximation
The cost of constructing and storing the tensors T and A from Section 2.1 grows exponentially with d.
To mitigate this growth, we replace T by a low-rank approximation T .The following lemma is a direct generalization of [28,Lemma 2.2] bounding the error introduced by such an approximation.
Lemma 1.Let f be defined as in (3). ) , consider the polynomial with where ∥•∥ ∞ denotes the uniform norm for functions and the maximum norm for tensors.
In the following, we approximate the tensor T by a tensor T in TT format [68], the discrete analogue of (2).The entries of T take the form with the so-called TT cores Given T in the TT format (6), the tensor Â defined in Lemma 1 can also be expressed in TT format with the TT cores A (ℓ) := G (ℓ) × 2 F (ℓ) .In turn, the function f defined in(5) can be expressed in the FTT format (2), with the corresponding univariate functions given by

Extended functional tensor train (EFTT) approximation
As discussed in the introduction, we can potentially compress FTT approximations further by essentially replacing the TT core G (ℓ) in ( 6) by H (ℓ) × 2 U (ℓ) , where To compute such a compressed approximation without explicitly forming and then compressing a FTT approximation, we propose to compute an approximation of the evaluation tensor (4) in the extended TT format [79].In this work, we obtain such an approximation by first approximating the evaluation tensor T from (4) in the Tucker format [94]: where )×r ℓ are called factor matrices.The multilinear rank of T is given by r 1 , . . ., r d (more precisely, the smallest values for these integers such that T admits ( 7)).
In a second step, the core tensor C is approximated in TT format: with the TT cores H (ℓ) ∈ R R ℓ−1 ×r ℓ ×R ℓ .Inserted into (7), this yields an approximation in extended TT format: This only requires O(drR 2 + dnr) storage, where R = max R ℓ , r = max r ℓ , n = max n ℓ , which compares favorably with the O(dnR 2 ) storage needed by the TT approximation (6) alone, especially when when r ≪ n.
Figure 1: Tensor network representation [66] of the coefficient tensor Â in (10) corresponding to EFTT format (11).The colored boxes mark the subtensors corresponding to the approximation of the evaluation tensor T as in (5) and the approximation of the core tensor Ĉ as in (8).

Approximation algorithm
In the following, we develop a novel algorithm for computing approximations in the EFTT format (11).Our algorithm obtains the factor matrices U (1) , . . ., U (d) in the Tucker approximation (7) from fibers of the evaluation tensor T via a variant of column subset selection [19,24].Following [28], the core tensor C is obtained as a subtensor of T by applying discrete empirical interpolation.Thus, there is no need to form C explicitly and we apply a variant of the TT-cross algorithm [67] to compute the TT cores H (1) , . . ., H (d) for the TT approximation (8) from only some entries of C.

Factor matrices
Fiber approximations.When fixing all but the ℓth index of T , one obtains a vector called a mode-ℓ fiber of X .The mode-ℓ matricization of T is a matrix that collects all mode-ℓ fibers [55].Applying adaptive cross approximation (ACA) [10], a popular low-rank approximation technique that accesses a matrix via its entries, to T {ℓ} yields index sets Îℓ , Ĵℓ that determine an approximation of the form Note that the matrix T {ℓ} (:, Ĵℓ ) contains mode-ℓ fibers of T .If the approximation error of ( 12) is small, we choose Û (ℓ) = T {ℓ} (:, Ĵℓ ) as an approximate basis of mode-ℓ fibers.
ACA with full pivoting determines the indices in Îℓ , Ĵℓ subsequently by choosing the entry of largest absolute value in the current approximation error.This corresponds to greedily maximizing the volume of T {ℓ} ( Îℓ , Ĵℓ ) [35] and comes with theoretical guarantees [20].However, it is impossible to apply full pivoting in our context because it requires the evaluation of all entries.Even partial pivoting [11] is impractical as it requires the evaluation of a (very long) row in T {ℓ} in every step.Inspired by the success of randomization in numerical linear algebra [59,97], we propose to determine the next indices based on sampling.This leads to Algorithm 1, which samples a fixed number s of random entries and picks the entry of largest absolute value.These random entries are also used for stopping the algorithm and determine the approximation rank r ℓ adaptively.Note that the update in line 10 is only performed implicitly and the subtraction is carried out each time an entry of A is evaluated.Overall, Algorithm 1 applied to T {ℓ} determines the index sets Îℓ , Ĵℓ using only O(r 3 ℓ + sr 2 ℓ ) entries of T .Afterwards, we explicitly compute Û (ℓ) = T {ℓ} (:, Ĵℓ ) by evaluating n ℓ r ℓ entries of T .We perform the described procedure for every ℓ = 1, . . ., d to determine Û (1) , . . ., Û (d) ; see also lines 4 and 5 of Algorithm 3 below.Construct S ⊂ {1, . . ., n} × {1, . . ., m} by uniformly sampling s index pairs (i, j).

Algorithm 1 Adaptive cross approximation with randomized pivoting
Tucker approximation To arrive at a Tucker approximation (7) we need to project the fibers of T onto the spans of )×r ℓ has orthonormal columns.The orthogonal projections used in the higher order singular value decomposition (HOSVD) [21], i.e. multiplying T with Q (ℓ) (Q (ℓ) ) T in each mode, would require to explicitly construct the whole tensor T .To avoid this, we follow the ideas of [28,82].Given a set I ℓ containing r ℓ indices selected from {0, . . ., n ℓ }, we let Φ I ℓ contain the corresponding columns of the (n ℓ + 1) × (n ℓ + 1) identity matrix, such that defines an oblique projection onto the span of Û (ℓ) .Applying this projections to each mode of T yields the Tucker approximation with factor matrices ) −1 and core tensor C = T (I 1 , . . ., I d ).Note that the subtensor C of T will not be explicitly formed; we only need to evaluate some its entries when computing the TT approximation of C in Section 3.2.We want to remark that the core tensor obtained using the HOSVD is not a subtensor of T .
In [28], it is shown how the error of the approximation ( 13) is linked to where ∥•∥ 2 , ∥•∥ F denote the spectral and Frobenius norms, respectively.The term ) T ∥ F is small when all mode-ℓ fibers of T are well approximated by the span of Û (ℓ) .The term ∥(Φ T I ℓ Q (ℓ) ) −1 ∥ 2 depends on the choice of the index set I ℓ .We apply discrete empirical interpolation [16] as formalized in Algorithm 2 to Q (ℓ) to find a suitable index set I ℓ .The whole process of computing U (ℓ)  and C from Û (ℓ) is summarized in Algorithm 3. Alternative randomized algorithms to compute Tucker approximations have been proposed in [63,75].Ĵℓ ← index set J returned by Algorithm 1 applied to T {ℓ} with tolerance ε and s samples.

5:
7: For simplicity, we have assumed that the values of n 1 , . . ., n d , determining the size of the expansion (3) and of T are given as input.In practice, these values are usually not provided.To ensure that the approximation error of the final approximation (11) is small, we use a heuristic to determine n 1 , . . ., n d adaptively.We initially set n 1 = • • • = n d = 16.After computing Û (ℓ) in line 5 of Algorithm 3, we apply the Chebfun chopping heuristic proposed by Aurentz and Trefethen [5] to Û (ℓ) .This heuristic decides based on the decay of the coefficients in the columns of the matrix F (ℓ) Û (ℓ) and the tolerance ε whether n ℓ is sufficiently large.If the chopping heuristic indicates that n ℓ is not sufficiently large then one sets n ℓ ← 2n ℓ + 1, updates T and repeats the procedure from line 4 to compute a new index set Ĵℓ .

TT cores
It remains to compute a TT approximation of the core tensor (8) to obtain the desired extended TT approximation (9).For this purpose, we use the rank-adaptive greedy restricted cross interpolation algorithm [78] as implemented in the TT toolbox 1 .This algorithm requires the evaluation of only O(drR 2 ) entries of C where R = max R ℓ , r = max r ℓ .
A practical algorithm is obtained by updating the index sets sequentially for ℓ = 1, . . ., d − 1 until we obtain a good approximation.For this purpose so-called DMRG supercores [77] are formed.These are defined as subtensors C(I ≤ℓ−1 , :, :, I >ℓ+1 ) which are reshaped into tensors in R R ℓ−1 ×r ℓ ×r ℓ+1 ×R ℓ+1 .The rank-adaptive algorithm proposed by Savostyanov [78, Algorithm 2] computes cross approximations of supercores reshaped into matrices in R (R ℓ−1 r ℓ )×(r ℓ+1 R ℓ+1 ) .After reshaping back into tensors, this leads to an approximation of the form where the index sets J ≤ℓ , J ≥ℓ are constructed such that I ≤ℓ ⊂ J ≤ℓ and I >ℓ ⊂ J >ℓ .This is achieved by initializing Algorithm 1 with the index sets I ≤ℓ , I >ℓ instead of empty sets.The sampling of the index tuples in Line 6 is modified to ensure the indices I ≤ℓ−1 and I >ℓ+1 are nested in J ≤ℓ and J >ℓ respectively.Moreover, Algorithm 1, is stopped after the first iteration, i.e. the index sets J ≤ℓ , J >ℓ contain at most one element more than I ≤ℓ , I >ℓ .We can now enhance the approximation ( 14) by replacing the index sets I ≤ℓ , I >ℓ by J ≤ℓ , J >ℓ .Note that this step increases R ℓ when I ≤ℓ ̸ = J ≤ℓ .The index update is repeatedly performed for ℓ = 1, . . ., d − 1 and stopped once the error of the approximation ( 14) is sufficiently small at sample points.We formalize this procedure in Algorithm 4.

EFTT approximation algorithm
We formalize the overall procedure for computing approximations in the EFTT format (11) in Algorithm 5. Note that r 1 , . . ., r d and R 1 , . . ., R d−1 are determined adaptively in lines 5 and 6 of Algorithm 5. Let R = max R ℓ , r = max r ℓ , n = max n ℓ .The total number of evaluations of f in Algorithm 5 is O(dr 3 + dnr + dsr 2 +drR 2 ).Applying Algorithm 4 directly to T requires O(dnR 2 ) evaluations and yields an approximation in FTT format (6).Note that Algorithm 5 may also be of interest in other situations, not necessarily related to functions.Lines 5-7 are of relevance whenever the entries can be (cheaply) evaluated but explicit storage of the tensor as a whole is too expensive.

Algorithm 5 EFTT approximation
number of samples s 2: Output: TT cores H (ℓ) and procedures for evaluating univariate functions u (ℓ) j defining an approximation of f in the EFTT format (11) 4: Define a procedure to evaluate entries of the evaluation tensor T defined in (4).5: Apply Algorithm 3 with tolerance ε and s samples to T to determine Û (ℓ) .Simultaneously, we update n 1 , . . ., n d as well as T as described in Remark 1. 6: Compute the tensors H (ℓ) by applying Algorithm 4 with tolerance ε to C. 7: Create procedures to evaluate univariate functions u (ℓ)

Numerical experiments
In this section, we present numerical experiments2 to study the performance of Algorithm 5. Unless mentioned otherwise the sample size in Algorithm 1 is set to s = min{n/2, 50}, where n = ((n 1 + 1)(n 2 + 1) • • • (n d + 1)) 1/d , and all tolerances are set to 10 −10 .The approximations error is measured via a Monte-Carlo estimation of the relative L 2 -error using 10 000 samples.

Comparison to a direct TT approximation
We first compare Algorithm 5 yielding an approximation in the EFTT format (11) to a direct TT approximation of the evaluation tensor T ∈ R 100×•••×100 (4) using Algorithm 4 as in [13].The latter approach yields an approximation in the FTT format (2).Since a direct TT approximation is not adaptive with respect to the polynomial degree, we fix the polynomial degree throughout this comparison, i.e., we replace line 3 by n 1 = • • • = n d = 100 and do not update adaptively in line 5 of Algorithm 5.
Benchmark functions.For the set of benchmark functions defined in Appendix A, we on average reduce the number of function evaluations required by 30.6% and the required storage by 41.6% when using Algorithm 5 compared to the direct TT approximation of the coefficient tensor.For the Ackley function the reduction is 88.8% in terms of function evaluations and 93% in terms of storage.For the Borehole function, we have r ℓ = R 2 ℓ , i.e. we can not compress the FTT approximation further.In this case, the direct TT approximation is slightly more efficient.At the same time, this demonstrates that all other test functions, taken from a wide range of applications, can be stored more efficiently using our approach.Details, including insights into the reliability of the algorithms, can be found in Table 5 in the appendix.Algorithm 5 requires fewer function evaluation compared to a direct TT-cross approximation.However, in very large dimensions (d > 300) the direct TT-cross approach leads to slightly more accurate approximations.A possible explanation for this observation might be that, for Algorithm 5, the error of d projections needs to be taken into account in addition to the TT-cross approximation error.Since the errors introduced by each projection might affect the overall approximation error multiplicatively (see [28,Theorem 2]), this could lead to an approximation error growing exponentially in very large dimensions.Moreover, approximating the smaller tensor C instead of T directly, might lead to slightly worse conditioned inverse matrices in (14).[13] and display the geometric mean of the error and the arithmetic mean for the number of function evaluations.

Comparison to the FTT approximation algorithm
In this section, we compare the performance of our proposed Algorithm 5 based on the EFTT format (11) to the approximation algorithm proposed in [36,38] and implemented in the c3py package 3 .The c3py algorithm uses a continuous variant of the TT-cross algorithm to compute approximations in the FTT format (2).The c3py package is based on the Legendre polynomials P 0 (x) = 1, P 1 (x) = x, (k + 2)P k+2 (x) = (2k + 3)xP k+1 (x) − (k + 1)P k (x) for k = 0, . . ., m − 2 [54] with degree at most m ∈ N instead of Chebyshev polynomials.This leads to multivariate polynomial approximations with degree at most (m 1 , . . ., m d ) of the form where ) .The package essentially computes B in TT format (6).The c3py package obtains the m coefficients for approximating a univariate function using Clenshaw-Curtis quadrature [18] in n > m points to approximately compute L 2 -projections onto the Legendre basis functions [36].To obtain a fair comparison [14], we modify our approximation algorithm slightly to also include such a Clenshaw-Curtis quadrature.Given the evaluation tensor d) , where E (ℓ) ∈ R (n ℓ +1)×(m ℓ +1) is defined as . . .
2 ) . . .2n ℓ +1 2 w (ℓ) , where the nonnegative weights w i ℓ are defined for i ℓ = 0, . . ., m ℓ , ℓ = 1, . . ., d as in [96]: Note that the matrix E (ℓ) encodes Clenshaw-Curtis quadrature in n ℓ + 1 nodes to (approximately) compute the Legendre coefficients via L 2 -projections as in c3py.Using analogous constructions as in Section 2, we can transform an extended TT approximation of the evaluation tensor (9) into a EFTT approximation with univariate functions represented in terms of linear combinations of Legendre polynomials.This can also be seen as storing B in (15) in the extended TT format (9).
In the following experiments we set n ℓ = 2m ℓ for ℓ = 1, . . ., d to ensure accurate quadrature.To determine the polynomial degrees m ℓ adaptively (see Remark 1), we follow the fiber adaptation strategy of c3py (see [36, Section 3.6.1]):We progressively increase the degrees until four sequential Legendre coefficients are smaller than the tolerance 10 −10 or until the maximum of m ℓ = 105 has been reached.In the c3py algorithm, we set all tolerances to 10 −10 .
Benchmark functions.We apply both the c3py algorithm from [38] and our novel Algorithm 5 to approximate the benchmark functions defined in Appendix A. Our numerical results displayed in Table 6 in the appendix demonstrate that our novel algorithm typically requires fewer function evaluations and still achieves the same accuracy compared to c3py.For the Ackley function, our approach reduces the number of required function evaluations by more than 96% and the storage by more than 90%.At the same time, our approach is slightly more accurate.

Integration of sin function
In the following, we repeat the experiment from [36, Figure 3-6].The function f (x 1 , . . ., x d ) = sin(x 1 + x 2 + • • • + x d ) can be represented in FTT format with TT ranks (2, 2, . . ., 2) and its integral over the domain [0, 1] d is known analytically [36].In Figure 3, we plot the error of the integral of the approximations computed via Algorithm 5 and the c3py algorithm.The figure shows, that our novel approach achieves a similar level of accuracy.The difference in the number of function evaluations is rather small, since the multilinear ranks and TT ranks of the function are small.Application: Uncertainty quantification.A classical application of multivariate function approximation is the computation of surrogates for uncertainty quantification [84,98].For the approximation of the quantity of interest mapping defined in Appendix C, we find in Table 1

Conclusion and outlook
For a broad range of functions arising in various applications, including mechanical engineering and uncertainty quantification, our novel algorithm for computing approximations in the EFTT format outperforms existing algorithms based on the FTT format in terms of function evaluations.In this paper, we have not discussed how to compute numerically with functions once they are compressed in EFTT format.For arithmetic operations, one would need to extend the ideas from [86,89] to EFTT.When a function is defined implicitly as the solution of a partial differential equation, extensions of [83,87] based on ideas from [57,70,80] x
The parameters w i and c i are drawn uniformly from [0, 1], where w i act as a shift for the functions while c i determines the approximation difficulty of the functions.We normalized c i such that where the scaling constants h and b are defined for each function in Table 4.Note that f 1 can be represented in FTT format (2) with max R ℓ = 2.The function f 3 is separable.For f 2 , we are not aware of any analytic FTT representation.4: Scaling constants for the Genz functions as in [13].

C Parametric PDE problem
In the following, we recall the example from [58,Section 4] where u(x, p) denotes the solution of the PDE (16) for the given parameter p ∈ [−1, 1] d .For each value of p, we solve the resulting PDE using a discretization based on linear finite elements.Table 5: As described in Section 4.1, we apply Algorithm 5 (EFTT) and an approximation of the evaluation tensor using Algorithm 4 (DirectTT) to approximate the test functions defined in Appendix A).For each function, we display the estimated L 2 error of the approximations, the number of function evaluations required to construct the approximation, the degrees of freedom in the approximation and the largest R ℓ and r ℓ of the resulting approximation.The numbers displayed are the mean (geometric mean for the error and arithmetic mean for all other quantities) after approximating each function 100 times.In addition, we display the standard deviation of the logarithm of the error σ(log(error)), as well as the standard deviation for the number of function evaluations (σ(#evals)) and the number of degrees of freedom (σ(#dofs)).(15).For each function, we display the estimated L 2 -error of the approximations, the number of function evaluations required to construct the approximation, the degrees of freedom in the approximation and the largest m ℓ , R ℓ and r ℓ of the resulting approximation.

Figure 2 :
Figure 2: We apply Algorithm 5 (EFTT) and an approximation of the evaluation tensor using Algorithm 4 (DirectTT) to approximate the Genz functions (see Appendix B) for dimension d ∈ {20, 50, 100, 200, 300, 400, 500}.Each row corresponds to one test function.Left: We plot the L 2 error for EFTT and DirectTT.Right: We plot the number of function evaluations required to compute the approximation for EFTT and DirectTT.Throughout this figure, we sample 30 different parameters for the Genz functions as in[13] and display the geometric mean of the error and the arithmetic mean for the number of function evaluations.

Figure 3 :
Figure 3: We apply Algorithm 5 (EFTT) and the algorithm in the c3py package to approximate f (x 1 , . . ., x d ) = sin(x 1 + x 2 + • • • + x d ) by functional low-rank approximation with Legendre polynomial basis functions (15) as in [36, Figure 3-6].Left: We plot the relative error of the integral of the approximations.Right: We plot the number of function evaluations required to compute the approximation.
d.If the maximal TT rank R = max R ℓ remains modest, this representation requires a storage of size O((d − 2)nR 2 + 2nR) with n = max n ℓ .This results in linear (instead of exponential) growth with respect d under the (strong) assumption that R remains constant as d increases.

Table 1 :
We apply Algorithm 5 (EFTT) and the algorithm in the c3py package to approximate the quantity of interest map Q : [−1, 1] d → R defined in Appendix C for d ∈ {4, 9, 16}.All tolerances are set to 10 −3 .The table displays the L 2 -error, the required number of evaluations (evals) and the number of degrees of freedom (dofs) for both approximations.
could be considered.

Table 6 :
dofs max ℓ m ℓ max ℓ R ℓ max ℓ r ℓ As described in Section 4.2, we apply Algorithm 5 (EFTT) and the algorithm in the c3py package to approximate the test functions defined in Appendix A by functional low-rank approximation with Legendre polynomial basis functions