Tensor rank bounds for point singularities in ℝ3

We analyze rates of approximation by quantized, tensor-structured representations of functions with isolated point singularities in ℝ3. We consider functions in countably normed Sobolev spaces with radial weights and analytic- or Gevrey-type control of weighted semi-norms. Several classes of boundary value and eigenvalue problems from science and engineering are discussed whose solutions belong to the countably normed spaces. It is shown that quantized, tensor-structured approximations of functions in these classes exhibit tensor ranks bounded polylogarithmically with respect to the accuracy ε ∈ (0,1) in the Sobolev space H1. We prove exponential convergence rates of three specific types of quantized tensor decompositions: quantized tensor train (QTT), transposed QTT and Tucker QTT. In addition, the bounds for the patchwise decompositions are uniform with respect to the position of the point singularity. An auxiliary result of independent interest is the proof of exponential convergence of hp-finite element approximations for Gevrey-regular functions with point singularities in the unit cube Q = (0,1)3. Numerical examples of function approximations and of Schrödinger-type eigenvalue problems illustrate the theoretical results.


Introduction
Recent years have seen the emergence of structured numerical linear algebra in scientific computing and data science. We mention only formatted matrix algebras, such as H-matrices (e.g., [33] and the references there) and tensor formats (e.g., [32,43,46,62,67] and the references there). To date, the impact of these methods was, first and foremost, on the corresponding scientific computing applications: being abstracted from fast multipole methods, formatted computational matrix algebras impact directly the numerical solution of elliptic and parabolic partial differential equations (see, e.g., [6,27,35]). Numerical tensor algebras, derived from quantum chemistry (e.g., [6,72] and the references there) have obvious applications in data science, where massive n-way data naturally arises and needs to be efficiently handled numerically. Furthermore, tensor-structured formats have, in recent years, been linked to deep neural networks (see [45,51] and the references there). We now comment on more specific developments in these areas which are directly related to the present paper, and the mathematical results obtained in it.
We are concerned with the approximation of functions with isolated point singularities using tensor-structured representations. In particular, we approximate, using quantized tensor decompositions, three-dimensional arrays of coefficients associated with the finite element projection of functions over trilinear Lagrange basis functions.
Quantization refers to the reshaping of an array of coefficients of size 2 × 2 × 2 into a multidimensional array of size 2×· · ·×2. The application of tensor decompositions (e.g., the Tensor-Train decomposition [65], which leads to the QTT-quantized tensor train decomposition, introduced in [42,64]) to such an array can lead to a reduction in complexity and number of parameters.
The number of parameters in a decomposition is related to the rank of the decomposition-i.e., the generalization of matrix rank to multi-dimensional arrays. Having a priori knowledge that a function of interest, e.g., the solution to a partial differential equation, can be approximated by a low-rank tensor decomposition, allows for the application of tensor-structured algorithms that avoid working with full 2 × 2 × 2 arrays of coefficients.
In particular, here we consider functions in weighted Sobolev spaces with radial weights and analytic-or Gevrey-type control of weighted semi-norms. Such functions arise in a variety of scientific applications: nonlinear Schrödinger equations (e.g., [11,13] and the references there), Hartree-Fock and density functional theory equations, continuum models of point defects [53], blowup solutions in evolution equations with critical nonlinearity (e.g., [70] and the references there) to name but a few.
The main result of the present paper is exponential convergence of tensorstructured approximations of point singularities in R 3 , i.e., they admit tensor ranks bounded polylogarithmically with respect to the accuracy ε ∈ (0, 1) of the approximation, measured in the Sobolev space H 1 .
An auxiliary result of independent interest is the exponential convergence of hpfinite element (FE) approximations for the class of functions considered. Due to the piecewise polynomial structure of hp-FE approximants, we can obtain their quantized representations with exact rank bounds that depend only on the dimensions of hp-spaces. This, in turn, leads to the desired rank bounds of the functions of interest.
One of the advantages of using quantized tensor decompositions-compared with the direct application of hp-FE approximations-is the relative ease of implementation. The adaptation of the number of parameters in the decomposition to the approximated function is based on well-known numerical linear algebra tools such as QR and SVD decompositions. Moreover, there exist open source codes with the implementation of basic linear algebra operations including solution of linear systems, which can be used independently of a particular application.
Note also that we do not need to know a priori the type and exact location of the singularity of the solution to solve PDEs in quantized tensor-structured formats. The nonlinear structure of the decomposition allows for an "automatic" adaptation of the tensor compressed representation to the regularity of the function. This is by contrast to hp methods, where mesh and polynomial degree refinements are programmed explicitly depending on the type of singularity. Furthermore, while the mesh of an hp space has to be constructed so that the refinement happens towards the singular point, this a priori knowledge is not necessary in the computation of quantized tensorstructured representation.

Tensor-structured function approximation
With the availability of efficient numerical realizations of tensor-structured numerical linear algebra, a new perspective has been opened towards computational function approximation. Here, one compresses arrays of function values in tensor formats; early work in this direction is [78], and [43] contains a bibliography with a large list of ensuing developments based on this idea. An (incomplete) list of references contains [37,42,44,61,66] where tensor rank bounds for specific functions have been obtained, both analytically and computationally, in the so-called quantized tensor train (QTT) format. QTT-formatted numerics for electron structure computations were presented in [41]. An analysis of approximation properties of tensor networks for classes of functions of finite differentiability as expressed by membership in classical Sobolev and Besov spaces of finite order has recently been presented in [1][2][3].
Subsequently, and more directly related to the present work, rather than rank bounds for individual functions, tensor rank bounds for solution classes of elliptic PDEs in one and two spatial dimensions were obtained in [36,38,40,58]. In [40], in particular, it was proved first that functions in countably normed, analytic function classes in polygons D ⊂ R 2 admit QTT-structured tensor approximations with tensor ranks bounded polylogarithmically in terms of the approximation accuracy ε.
The key mathematical argument in the references cited above is based on analytic regularity results for solutions of elliptic PDEs in polytopal domains. Such regularity results, implying solutions belong to countably normed spaces, have been obtained in the past two decades for several broad classes of (boundary value and eigenvalue problems of) elliptic PDEs, in [4,55,59].

Contributions
In this work, we obtain rank bounds for the QTT approximation of functions in weighted analytic and Gevrey classes with point and edge singularities. Specifically, we prove polylogarithmic growth of QTT ranks with respect to the accuracy of the approximation. To obtain the bounds, we approximate functions from the weighted Gevrey class by continuous, piecewise polynomial functions on dyadic partitions which are geometrically refined towards the singular supports. These piecewise polynomial approximants are subsequently re-interpolated and compressed to the QTT format. The resulting rank bounds follow from the low-rank structure of piecewise polynomial functions, and the exponential convergence of the piecewise polynomial approximations.
The principal novel contributions of this paper are therefore the exponential convergence of hp approximation for weighted Gevrey classes, and the polylogarithmic rank bounds of QTT approximations on the these classes.
First, we analyze approximation rates of tensor-structured approximations of smooth functions with isolated point singularities. As compared to exponential convergence results for analytic functions with point singularities, we here establish exponential convergence of hp-finite element (FE) approximations on geometric meshes of axiparallel quadrilaterals resp. hexahedra analogous to [20] also for Gevrey-regular functions.
We then address tensor-formatted approximations. Generalizing results also in two variables, in the present paper we extend the analysis in [40] to quantized, TTstructured function approximation of functions from countably weighted, Gevreytype classes. The corresponding results in three spatial variables are novel. They also extend the QTT rank bounds in [40] to Gevrey-d regular functions (see Remark 4). They also constitute a building block for the derivation of corresponding QTT rank bounds for edge and face singularities in three space dimensions, which we do not detail here. In particular, we prove in three physical variables for analytic and Gevrey functions with point singularities, for the classical tensor-format, asymptotic upper bounds on quantized tensor ranks at prescribed accuracy ε which are better than the corresponding bounds for the transposed TT format introduced in [40] (in two dimensions).
We show numerical results indicating the correctness of the presently obtained results, and also strongly suggesting that similar ranks are achieved in tensorformatted PDE solvers, provided the PDE solutions belong to the countably normed classes introduced in Section 1.3.2.

Problem formulation
The tensor-formatted function approximation considered in this paper aims at establishing tensor rank bounds for functions in certain classes of locally smooth functions that admit a point singularity. In this paper, we confine ourselves to the case that the function under consideration admits singular support consisting only of one isolated point (we therefore speak of "point singularities"). Naturally, functions whose singular support comprises of a finite number of well-separated points can equally be approximated in the tensor-formats discussed here, with the same tensor rank bounds, by a localization and superposition argument.
Weighted Sobolev spaces for functions with isolated point singularities have been introduced for the analysis of elliptic problems in polygonal domains, see [47], since they allow for the extension of classical elliptic regularity theory to domains with corners. For an overview of regularity results for elliptic boundary value problems in conical domains, in weighted Sobolev spaces, we refer to the monographs [28,48,49,60].
For elliptic boundary value problems in three space dimensions, weighted Sobolev spaces that accommodate isolated point singularities have also proven important in the mathematical regularity analysis of problems with singular potentials, such as electron structure calculations in quantum physics and quantum chemistry, see, e.g., [21][22][23].
When a function is regular in weighted Sobolev spaces-specifically, when analytic-type bounds can be derived on the norms of its derivatives-piecewise polynomial approximations can be constructed, for example by hp finite elements which converge exponentially (in terms of the number of parameters) [29,30,71,73,74]. This suggests the existence of an underlying low-rank structure in suitable tensor formats; for this reason, we are here interested in the derivation of rank bounds for functions that belong to weighted analytic-and Gevrey-type classes.
A theory of analytic regularity in weighted Sobolev spaces has been developed for several classes of important physical problems and we mention an incomplete list. Solutions to scalar elliptic problems with constant coefficients belong to analytictype weighted spaces [15,16], as do the flow and pressure obtained with the Stokes [31] and Navier-Stokes [59] equations in polygons. Furthermore, eigenfunctions to three-dimensional linear [55] and nonlinear [54] Schrödinger equations are weighted analytic. In quantum chemistry, the wave functions computed with the nonrelativistic Hartree-Fock models for electronic structure calculations are also analytic in weighted Sobolev spaces [56,Section 7.4], [12], with point singularities at the nuclei. We refer to Section 1.3.3 for some explicit examples in this sense. Other instances of the occurrence of point singularities in otherwise smooth solutions comprise general relativity (see, e.g., [10,79] and the references there) and solutions of parabolic evolution equations with critical nonlinearity (see, e.g., [70] and the references there). The results of the present work apply to all the problems cited above, whose solutions are weighted analytic, in the spaces that we detail in Section 1.3.2 below.
We consider the following setting for quantized, tensor train (TT)-formatted function approximation in Q = (0, 1) 3 , with one point singularity at the origin, where the functions belong to countably normed, weighted Sobolev spaces, where the weights are powers of r = |x|, the Euclidean distance of the point x ∈ Q from the origin.

Kondrat'ev-type weighted Sobolev spaces
For integer s ≥ 0, a real parameter γ ∈ R, and summability exponent 1 ≤ q < ∞, we introduce the homogeneous weighted Sobolev spaces and norm Our focus will be mostly on non-homogeneous weighted Sobolev spaces (remark the different weight exponent) In the following, we will always consider the case where q = 2, 0 < γ − 3/2 < 1, and s > γ − 3/2. Under those hypotheses, as shown in [14,Proposition 3.18], the above norm is equivalent to Non-homogeneous spaces allow for functions with non trivial Taylor expansion at the singularity and have been used, for this reason, in the analysis of problems in non smooth domains with Neumann boundary conditions and of elliptic problems with singular potentials. For a thorough analysis of the relationship between homogeneous and non-homogeneous spaces, we refer the reader to [49] and [14].

Gevrey and analytic function classes
We define the weighted Kondrat'ev-type class of functions of infinite regularity (Q) Furthermore, for constants C, A > 0 and d ≥ 1, we introduce the countably normed, homogeneous weighted Gevrey-type (analytic-type when d = 1) class The countably normed, non-homogeneous weighted classes J ∞,q γ (Q) are then defined as in the homogeneous case, while the non-homogeneous Gevrey/analytic classes are given by ; similarly we omit the summability exponent q when it equals 2 in the notation for the weighted Gevrey classes.

Model problems
We illustrate the scope of problems by listing several concrete boundary-value and eigenvalue problems whose solutions are known to belong to the weighted analytic classes K γ (Ω) and J γ (Ω). Although the focus here is on three-dimensional problems, we start by considering a polygon Ω ⊂ R 2 with n ≥ 3 straight sides and corners c i , i = 1, . . . , n. In this setting, the space K γ (Ω) contains the corner weight function r P = n i=1 |x − c i |, i.e., the seminorm (1) is replaced, for 1 ≤ q < ∞, by Then, given an analytic (in Ω) external force field f , the Stokes equations −ν u + ∇p = f in Ω, ∇ · u = 0 in Ω and the viscous, incompressible Navier-Stokes equations with homogeneous Dirichlet ("no-slip") boundary conditions have been shown in [31,59] to admit solutions in K γ (Ω) with γ > 3/2. Specifically, for the nonlinear boundary value problem (4) we require a "small data assumption" which is well-known to ensure uniqueness of Leray-Hopf solutions, see, e.g., [25,Chapter IV,Theorem 2.2]. See Remark 4 for further comments on the implication of the present work on two-dimensional problems.
In the three-dimensional setting, energy minimization problems in quantum physics/chemistry can be transformed into eigenvalue problems whose solutions are in the weighted analytic class (3). We consider here a set of isolated point singularities situated at n nuclei in positions R i ∈ R 3 , i = 1, . . . , n, and function spaces with weight function r such that r |x − R i | in the vicinity of each R i , and r 1 far from all singularities resp. all nuclei.
A first example is given by a nonlinear Schrödinger equation with polynomial nonlinearity. Consider a compact domain without boundary Ω (e.g., a periodic unit cell) and a potential V such that there exists β < 2 and a constant A V > 0 such that Then, the eigenfunction u corresponding to the smallest eigenvalue (i.e., the "ground state") of the nonlinear Schrödinger equation is in J γ (Ω) for some γ > 3/2, see [56,Section 7.3]. Note that (5) is the Euler-Lagrange equation of the minimization problem As a second example we consider the Hartree-Fock equation. Let V C be the potential of the Coulomb interaction exerted on electrons by nuclei with charge Z i assumed to be pointlike and situated at positions R i ∈ R 3 , i = 1, . . . , n, i.e., The Hartree-Fock model consists in finding the smallest N eigenvalues ε i and the corresponding Then, under some conditions on the potential V C so that the solution exists [52], the eigenfunctions are weighted analytic: see [56,Section 7.4]. Problem (6) is the Euler-Lagrange equation of the minimization problem (see [12,Section 9]) Remark 1 (Near-Singularity) While functions of the form (r ∈ R + , ω ∈ S 2 spherical coordinates) u a (r, ω) = (r 2 + a 2 ) β/2 v(ω), v analytic in S 2 (7) are, for real, nonzero values of a and for β > 0, formally (mathematically) smooth, their behavior approaches that of functions with point singularities at the origin when |a| 1. Specifically, if |a| ≤ a max , there exist positive constants C and A independent of a such that u a ∈ J γ (Q; C, A, 1) for γ < β + 3/2; hence, the bounds obtained in the present paper allow for the derivation of rank bounds for the quantized tensor-formatted approximations considered, which are uniform as the parameter a ↓ 0 for functions of the form (7).
The same remark applies to certain merging point singularities as arise, for example, in binary star or black hole models. Consider, e.g., two nuclei situated at locations R 1 = −εe 1 , R 2 = εe 1 in R 3 at distance 2ε for small ε > 0. Denoting by r i = |x − R i |, i = 1, 2, and r = |x|, we find v(x) = r 2 1 + r 2 2 = 2(r 2 + ε 2 ) i.e., once more a function of the above form with a = ε.

Structure of this paper
In Section 2, we review the definitions and notation of quantized, tensor-structured function approximation which are to be employed throughout the remainder of the article, extending the concepts of [66]. In Sections 2.2-2.7, in particular, we introduce the tensor train (TT), the quantized TT format (QTT), transposed quantized TT format (QT3) and the Tucker quantized TT format (TQTT), some of which allow to prove better rank bounds on functions with point singularities.
Section 3 introduces tools from numerical analysis which we require in the arguments for the TT rank bounds for function approximation. Section 3.1 introduces in particular the notion of "uniform background mesh" (never directly accessed in the QTT formats) which is the basis for all quantized TT function representations. Section 3.2 recapitulates several notions and auxiliary results from the theory of socalled hp-approximation from [73,74,77]. Section 4 introduces a combined (quasi) interpolation projector, which was introduced in [40] (in two dimensions) and which is crucial in establishing the rank bounds. Section 5 then contains statements and proofs of the main results of the present paper: tensor rank bounds for generic functions in the various countably normed classes introduced in Section 1.3 above. These bounds are obtained for functions with the singularity at a corner of the domain; they are extended to the case of an internal singular point (and to a patchwise formulation that allows for more complex domains) in Appendix B.
Section 6 presents detailed numerical experiments which exhibit actual TT rank bounds in the various formats for model singular functions in three space dimensions. The Section 7 provides a brief summary of the main results, and possible further research directions. Appendix A contains (novel) auxiliary results on exponential rates of convergence of hp-approximations for Gevrey-regular functions in R 3 with point singularities, generalizing [20] to axiparallel geometric meshes of hexahedra with 1-irregular edges and faces.

Tensor-structured representations
The mathematical issue in tensor-formatted function approximation consists in finding a compressed representation/approximation of three-way tensors ∈ N. All techniques that we examine are based on the Quantized Tensor Train (QTT) representation, see, e.g., [38,40,42,64,66] and the references there.
In particular, we will analyze three tensor compressed representations, that we call here (classic) QTT, transposed QTT (QT3), and Tucker QTT (TQTT) representation, respectively. The difference between these schemes lies in the arrangement of the three physical dimensions of the tensor A in the corresponding TT format in the following, after a brief introduction of QTT representations, we detail the three formats mentioned.

Notation
Throughout, we adopt the following notation, from [40]. Given n ∈ N indices i 1 , . . . , i n such that i j ∈ {0, . . . , k j − 1} for all j = 1, . . . , n, we write In what follows, the term tensor will generically denote a multi-dimensional array. Furthermore, for an axiparallel d-dimensional (d ≤ 3) subset K ∈ Q, the space Q p (K) is the tensor product space of d-variate polynomials in K of maximum polynomial degree p in each variable. Furthermore, we will indicate by a colon ":" a whole slice of a tensor. For example, given a four-dimensional tensor A ∈ R n 1 ×n 2 ×n 3 ×n 4 with entries a i,j,k,l , we will write A i,:,:,l = {a i,j,k,l } j =1,...,n 2 ,k=1,...,n 3 ∈ R n 2 ×n 3 .

Tensor train (TT) format
Tensor Trains (TT) [65], also known as Matrix Product States (MPS) in the computational physics community [72], provide an efficient way to represent highdimensional tensors, provided these tensors have an underlying low-rank structure. Let d 1, and consider the d-dimensional tensor The Tensor Train representation of the d-variate tensor B in (8) is given in terms of the core tensors 1 where r k ∈ N (with the restriction r 0 = r d = 1) and such that [65, Eq. (1. 2)] Suppose for ease of presentation that r i = r and n i = n for all i. Then, the TT representation (10) has N dof = O(dnr 2 ) parameters. The TT format is therefore an efficient decomposition if a d-way tensor can be written as a tensor train with low ranks r i . Due to the equality in (10) the representation we have introduced is an exact TT representation; in practice, a matrix may not admit an exact low-rank TT representation, but a low-rank approximation could instead be available.

Rank bound analysis of TT representations
To examine the issue of low-rank approximation of high-dimensional tensors, we require the concept of unfolding matrices ("unfoldings" for short), used to derive rank bounds on the TT representation of a tensor.
Definition 1 (Unfolding matrix) Let d ∈ N and n i ∈ N for i = 1, . . . , d. Given a tensor B ∈ R n 1 ×···×n d , we define for all q = 1, . . . , d − 1 its unfolding matrices B (q) ∈ R n 1 ···n q ×n q+1 ···n d as for all i k = 1, . . . , n k and k = 1, . . . , d i.e., the matrix with row index given by the concatenation of the first q indices, and column index given by the concatenation of the remaining ones.
In the case that the unfolding matrices of a tensor can be approximated by lowrank matrices, then a low-rank TT approximation exists. This is made precise in the following result.

Proposition 1 [65, Theorem 2.2]
Let B ∈ R n 1 ×···×n d such that its unfolding matrices B (q) can be decomposed as There exists a tensor C with TT representation (10) and TT ranks r q such that notation, the TT decomposition of B can be written as which is also used in [65, Eq. (1. 3)]. For clarity of presentation, we use the representation (10), which is more compact than (9) (and equivalent to it). We also do not distinguish between the mappings U k and the three-dimensional arrays U k .
The above theorem includes as a sub-case the rank bound of exact TT representation, by affirming the existence of an exact TT rank r q representation of a tensor with unfolding matrix rank bounded by r q , q = 1, . . . , d − 1.

Quantized tensor train (QTT) format in one physical dimension
We introduce QTT representations in the simplified setting of QTT approximation of vectors for ∈ N. Generalizations to the multi-dimensional case will be the subject of the next sections.
The QTT decomposition introduced in [42,64] extends the use of the TT approximation to the case of low-dimensional tensors with a large number of elements. To do so, the low-dimensional tensor is reshaped into a high-dimensional one, which is subsequently TT-(re)approximated. Applied to the vector in (11), algorithmically this is achieved by reshaping it into the -dimensional tensor v such that The tensor v can then be represented in TT form. We formalize this representation in the following definition.
Definition 2 (Univariate QTT decomposition) Given ∈ N and a vector v ∈ R 2 , v admits a QTT representation with QTT ranks r 0 , . . . , r and QTT cores As before, the tensor cores can also be interpreted as three-way arrays in R r i−1 ×2×r i .

Classic QTT format in three physical space dimensions
The "classic QTT" format is the straightforward generalization of the univariate QTT format in Definition 2 to the multivariate case.
In this way, a three-dimensional tensor A ∈ R 2 ×2 ×2 is reshaped into the tensor for all i n , j n , k n ∈ {0, 1}, which is subsequently TT-decomposed.

Transposed order QTT format in three physical space dimensions
In the transposed order QTT format (referred to as "QT3" format) introduced first in [40], after reshaping the tensor A as in (12), the indices from the different (three) physical dimensions are regrouped together, resulting in a tensor for all i n , j n , k n ∈ {0, 1}. The tensor A qt3 is subsequently TT-decomposed, as specified in the next definition. . Each node represents a tensor with edges in the network indicating indices. An edge connecting two nodes is a contracted index (corresponding to tensor multiplication). This can be seen comparing the networks with equations (13), (17), and (20) Definition 4 (Transposed order QTT) Let ∈ N and let A ∈ R 2 ×2 ×2 . The tensor A admits a transposed QTT decomposition with tensor ranks r 0 , . . . , r and cores for all i n , j n , k n ∈ {0, 1}, and where U n : {0, . . . , 7} → R r n−1 ×r n , for n = 1, . . . , .
We have the restriction on the ranks r 0 = r = 1.
A representation of the transposed order QTT decomposition in tensor network format is given in Fig. 1b.

Tucker QTT
The Tucker QTT (TQTT) decomposition is a combination of the Tucker and the QTT decompositions, first considered in [18] where R 1 , R 2 , R 3 ∈ N are the Tucker ranks, the tensor G ∈ R R 1 ×R 2 ×R 3 is the Tucker core and the Tucker factors U, V , W can be considered as matrices U ∈ R 2 ×R 1 , V ∈ R 2 ×R 2 , W ∈ R 2 ×R 3 . In the TQTT decomposition, the factor matrices U, V , W are given by QTT decompositions, where, e.g., for U , only one of the QTT cores depends on the corresponding column number β 1 : We denote the QTT ranks of U, V , W as {r 0 , r 1 , . . . , r }, {s 0 , s 1 , . . . , s } and {t 0 , t 1 , . . . , t } with the constraints r 0 = R 1 , s 0 = R 2 , t 0 = R 3 and r = s = t = 1.

Degrees of freedom
Supposing for ease of notation that r n = r qtt for the classic QTT representation, r n = r qt3 for the transposed one, and r n = s n = t n = r tqtt and R 1 = R 2 = R 3 = R for Tucker QTT, the number N dof of parameters in the QTT representations is bounded as

Functional setting
Our analysis will require the introduction of two different meshes and of two respective finite element spaces in the cube Q. The first one is a uniform tensor product mesh with distance between nodes given by h = 2 − . This mesh contains 2 nodes in every physical direction; given a function f defined over Q, the point values of f at the mesh points can be grouped in a three-dimensional tensor of dimension 2 × 2 × 2 , which can be QTT-approximated in the formats introduced in the previous section. Note that, in practice, one does not need to compute the values of the function at all 2 3 mesh points, see, e.g., [63], as this would undermine the efficiency of tensor compressed methods. For this reason, the background mesh is also referred to as "virtual" mesh in the literature, see for example [38][39][40]68]. Furthermore, tensor-formatted closed forms of some discrete differential operators exist, see, e.g., [37,42,44]. This can be used to discretize certain partial differential equations in quantized tensor format, as it will be shown in the sequel. The space of (tensor-formatted) functions on the uniform mesh is the space of Q 1 finite elements, i.e., the tensor product of one-dimensional Lagrange functions associated with mesh nodes.
The second finite dimensional space we introduce is the auxiliary hp space. This space is introduced here only for proving tensor rank bounds of the QTT-structured approximation. It is never accessed during numerical computation in the tensor formats. The hp space is, in particular, an H 1 -conforming finite element space, on a mesh with elements geometrically refined towards the origin. The polynomial degree of functions in the hp space is, instead, increasing polynomially with the number of geometric mesh layers. This is made more precise in Section 3.2 below. The role of the auxiliary hp finite element approximation is to provide an exponentially convergent, continuous and piecewise polynomial approximation on a bisection geometric partition which is compatible with the background mesh in the Q 1 -approximation for generic functions in the weighted Sobolev space J γ (Q).
A function in v ∈ J γ (Ω) can then be approximated-with exponential accuracy-by its projection v hp into the hp space. By re-interpolating v hp on the background mesh and QTT-compressing the resulting tensor, we establish existence of quantized, tensor-structured approximations with polylogarithmic bounds on the QTT ranks and the number of QTT parameters. The quasi-interpolation operator from v to its representation on the background mesh is introduced in Section 3.3.
For simplicity, we will consider here functions that have zero trace on the part of the boundary not abutting at the origin, i.e., on We denote by H 1 Γ (Q) the subspace of H 1 (Q) functions with zero trace on Γ . We then fix γ ∈ R such that γ − 3/2 ∈ (0, 1), two constants C X , A X > 0, and a regularity exponent d ≥ 1 and denote by the weighted space of Gevrey-d-regular functions with zero trace on Γ that will be considered henceforth.

Low order background FE space X
We introduce the so-called "background" (sometimes also referred to as "virtual") FE space discussed above. Here, it will consist of the space of continuous, piecewise trilinear functions on a uniform mesh of axiparallel hexahedral elements of size 2 − (so-called Q 1 -FEM), which we now introduce.

Uniform background mesh T
In Q = (0, 1) 3 , we introduce the uniform mesh T with nodes

Background finite element space
For (i, j, k) ∈ {0, . . . , 2 − 1} 3 , we denote by ϕ i,j,k the locally trilinear, continuous nodal Lagrange functions which satisfy where δ im denotes the Kronecker delta symbol for indices i and m.
The space of continuous, locally trilinear Lagrange functions on the (background) Note that the basis functions ϕ i,j,k and the space X are algebraic tensor products of the corresponding univariate functions, resp. spaces. We remark that for every v ∈ X holds v | Γ = 0.
Remark 2 X contains functions that vanish on Γ . We limit ourselves to this case for simplicity of notation; the extension of our analysis to functions with nonzero trace on Γ involves additional technicalities. We refer to [40] for the two-dimensional case.

Lagrange interpolation operator I
We denote by I the Lagrange interpolation operator on the uniform tensor mesh T . I.e., I : C(Q) → X is defined as

Analysis and synthesis operators
For ∈ N, A : X → R 2 ×2 ×2 and S : R 2 ×2 ×2 → X are the analysis and synthesis operators, such that

Auxiliary hp space
We obtain the QTT rank bounds on TT-formatted approximations by comparison with hp-approximations. To this end, we introduce the hp-FE spaces. We start with 1-irregular meshes of axiparallel hexahedra with geometric refinement towards the singularity of the function of interest ("geometric meshes" for short). Element K 000,0 has one vertex coinciding with the origin. We collect all elements at the same refinement level in mesh layers We also introduce the one-and two-dimensional versions of the geometric mesh as  We remark that for q ∈ N, k = 1, . . . , 2 q − 1, j = 1, . . . , q − 1, and for all integer

hp space
The hp space is formally introduced as for all n ∈ N , j = 1, . . . , and n = 000, j = 0}. (25) Note that, as a consequence of the existence of a continuous hp approximation to functions in J ∞ γ (Q) proved in Appendix A and in [73], the space X ,p hp is well defined by (25).

hp approximation
We provide a brief presentation of (novel) hp-interpolation error bounds which are exponential in the number of degree of freedom for functions in the Gevrey-type We consider here axiparallel, geometric partitions of Q = (0, 1) 3 into hexahedral elements; this entails, of course, irregular nodes and faces so that hp-interpolants are to be constructed in a two-stage process: first, an elementwise hp-(quasi)interpolant with analytic error bounds and second, polynomial face jump liftings which preserve the analytic bounds. We refer to the Appendix and to, e.g., [73] for details on this.
In the analytic case, i.e., when d = 1 such exponential error bounds are wellknown (e.g., [73]). However, for d > 1, these bounds are novel; for regular geometric meshes of tetrahedra, corresponding bounds have recently been established in [20].
We introduce in Appendix A the projector Π ,p and dim(X ,p hp ) p 3 3d+1 , see Proposition 3 in Appendix A.

Quasi interpolation operator P
We recall that x 1 x 2 x 3 = 0} and for γ > 3/2, positive constants C X and A X , and Gevrey exponent d ≥ 1. We also fix p d such that (26) holds and define the quasi-interpolation operator P : X → X as

Quasi interpolation error
We give here (specifically, in Proposition 2) an estimate on the error introduced by the quasi-interpolation operator P defined in Section 3.3. We start by estimating, in the following lemma, the error introduced by interpolating the hp projection of a function in X.
Proof There holds where the mesh layers L j are defined in (23). The quantity on the right-hand side of this equation is an upper bound for the error of interpolation over a uniform mesh of axiparallel cubes of edge length 2 − . The axiparallel cubes K n,j are affine equivalent to the reference element K = (−1, 1) 3 . Furthermore, since u ∈ X, by Remark 8 in the Appendix there holds Π ,p Hence there exists C > 0 such that, for all (n, j ) ∈ N × {1, . . . , } ∪ (000, 0) and for all By the polynomial inverse inequality where v is a polynomial of degree p and h K is the diameter of K (see, e.g., [24,77]), recalling that an element K n,j ∈ L j is an axiparallel cube with diameter h j 2 − +j ) and using a triangle inequality, there exists a constant C > 0 independent of and of p such that Since γ > 1, there exists a uniform constant C such that, on each K n,j , From this last inequality and (2) Combining (28) to (29), there exists a constant C > 0 such that for all holds Then, by (26) and since p d , Absorbing the terms algebraic in into the exponential by a change of constant concludes the proof.
where P u ∈ X is defined in (27) and there exists C > 0 independent of ε such that Proof By a triangle inequality, Lemma 1, and Proposition 3 (recalled above in (26)), where C and b are independent of . The choice = b −1 log C ε and adjusting the value of C concludes the proof.

QTT-formatted approximation of u ∈ J γ (Q)
We now state and prove our main results. For Gevrey-regular functions in Q = (0, 1) 3 with point singularity at the origin, and for each of the three tensor formats (QTT, QT3, TQTT), we prove bounds on the ranks which are sufficient to achieve a prescribed approximation accuracy ε ∈ (0, 1) in the norm H 1 (Q). This is the relevant norm for linear, second order, elliptic PDEs.
Proof We consider the unfolding matrices V (q) of T qtt (A (v qtt )), with T qtt defined in (15) and A in (22). We first consider the case q ∈ {1, . . . , − 1}. In this case, we denote its left and right endpoints as y K 0 and y K 1 , so that (y K 0 , y K 1 ) = K. On S 1 , we introduce a geometric mesh , and the univariate discontinuous FE space and v is right continuous at the nodes of G q 1d .

Consider the affine transformations
Then, x ξ,η 1 η 2 ,ζ ∈ ψ ξ,η 1 (Sq 2 ). Moreover, for each ξ, η 1 there exists a piecewise for all η 2 = 0, . . . , 2 −q − 1 and ζ = 0, . . . 2 − 1, see Fig. 4. Since dim(Xq S 2 ) ≤ C p 2 , we obtain, reasoning as before, It remains to consider q with 2 ≤ q < 3 . We sketch their treatment, which follows the same line of reasoning as in the preceding cases. Every row of the unfolding matrix V (q) contains the evaluation of v qtt on 2 3 −q equispaced points belonging to a line parallel to the z-axis. Hence, there exists a space of piecewise polynomials with less than C(3 − q)p degrees of freedom such that each row of V (q) can be written as linear combination of elements of the space, thus implying the existence of a constant C > 0 such that The proof is concluded by remarking that choosing | log ε| and using Proposition 2. Remark 3 In the proof of Lemma 2, we use a slightly different technique to treat the cases q ∈ {1, . . . , } and q ∈ { + 1, . . . , 3 }: in the first, we show that the columns of the unfolding matrices are point evaluations of piecewise polynomials, while in the latter we show a similar result for the rows of the unfolding matrices. It is possible to prove bounds on r 1 , . . . r by using the same procedure as the one used for q ≥ 2 + 1: one would introduce the three-dimensional slice and the space Then, using the affine transformation which implies by the same reasoning as used in the proof that for q ∈ {1, . . . , } This bound is weaker than the one obtained in the proof of Lemma 2 (and would impact the final bound on N dof ). The same effect would manifest itself if we changed the direction of the decomposition in the case q ∈ {2 + 1, . . . , 3 }.
Note that for ξ 1 = η 1 = ζ 1 = 0, the function v qt3 •ϕ ξ 1 ,η 1 ,ζ 1 is not a polynomial, which increases the dimension of the row space of U (q) by 1. Since dim(X S q ) ≤ Cp 3 , and using the same arguments as in the proof of Lemma 2 it can be concluded that

Fig. 5
Reference cube S q for transposed order QTT and representation of ϕ ξ 1 ,η 1 ,ζ 1 which gives, due to (21), a total number of degrees of freedom N dof ≤ C 6d+1 . The fact that | log ε| concludes the proof.

Rank bounds of TQTT approximations
In this section, we prove rank bounds for the TQTT approximation. We start by proving, in Lemma 4 below, rank bounds for the block QTT decomposition of collections of piecewise polynomial functions. Block QTT decompositions are precisely defined in the following.
We also need the definition, on the geometric mesh G 1d (see (24)), of the univariate hp-FE space Proof We provide a constructive proof with explicit formulas for the QTT cores. In the proof, the m-by-n matrix with zero entries will be denoted by O m×n , while the n-by-n identity matrix will be written as I n . Let Then q admits the exact, explicit QTT representation [26,42,66] with ranks r k = p + 1: Let now w α ∈ X We conclude that w α (x 0...01i k+1 ...i ) are polynomials sampled in equidistant points. By utilizing this fact, explicit formulas (31) for the polynomial parts and combining expressions for each of them, we obtain: and where we used the notation For fixed i k , k = 2, . . . , − 1, the matrices W k (i k ) are of size (p + s + 1) × (p + s + 1). We conclude that the ranks of the collection {W α } α are bounded from above by p + s + 1. This completes the proof.
Applying Lemma 4 to the Tucker factors U, V , W , we obtain their block QTT representation with tensor ranks {R, r QTT , . . . , r QTT , 1} bounded as To store Tucker factors in the block QTT format, we need to store N factors dof = 3 · 2 · (r QTT R + r 2 QTT ( − 2) + r QTT ) ≤ C 2d+3 parameters of the decomposition. Storing the Tucker core requires an additional N core dof = R 3 ≤ C 3d+3 parameters. This gives the following bound for the overall number of degrees of freedom in the TQTT representation Choosing | log ε| and using Proposition 2 completes the proof.

Exponential convergence of QTT approximations of u ∈ J γ (Q )
From Proposition 2 and Lemmas 2, 3, and 5, we obtain the following estimate for the QTT-Finite Element approximation of functions in J γ (Q). In the following theorem, we introduce a tag qtd ∈ {qtt, tqtt, qt3}, which generically denotes quantized tensor decomposition.
Theorem 1 Assume γ > 3/2, C u > 0, A u > 0, d ≥ 1, and 0 < ε 0 1. Furthermore, assume the function u belongs to the weighted Gevrey class u ∈ J γ (Q; C u , A u , γ, d) ∩ H 1 Γ (Q). Then, there exists a constant C > 0 such that, for all 0 < ε ≤ ε 0 , there exists ∈ N and v qtd ∈ X such that and v qtd admits a representation with Remark 4 (Rank bounds of QTT-formatted approximations of two-dimensional corner singularities) Using the same techniques for the two-dimensional case (which was already considered for the transposed order QTT in [40] in the analytic class, i.e., for d = 1) results in the bound κ = 2d + 3 for classic QTT, 4d + 1 for transposed order QTT.
In the case of two spatial variables, the Tucker QTT format is easily reduced to the classic QTT format for the index ordering i , . . . , i 1 , j 1 , . . . , j .

Numerical experiments
In this section, we support the obtained theoretical results with numerical experiments. First, in Section 6.1, we construct FE approximants to functions defined in Q = (0, 1) 3 with point singularities in three quantized tensor formats: QTT, transposed order QTT (QT3) and Tucker QTT (TQTT), see Fig. 1 for their tensor network representations. We note that for all formats under consideration, the numerically observed asymptotic behavior of rank versus error is better than that of theoretical estimates. In Section 6.2, we consider an elliptic eigenvalue problem with a singular potential -the Schrödinger equation for a hydrogen atom. We approximate the solution using the finite element method with a tensor of coefficients represented the TQTT format. The numerical results suggest that convergence rates of QTT-formatted approximations are slightly higher than those achieved by the hp-FEM.
We note that optimization algorithms, e.g., SVD, that we use to approximate tensors are based on optimization in L 2 -type norms. Nevertheless, in numerical experiments, we also calculate H 1 errors for function approximation problem (Section 6.1) as well as eigenvalue errors (Section 6.2) that are associated with the computation of derivatives. We observe that these errors behave in agreement with the theoretical bounds.

QTT-FE approximation of functions with point singularities
In this section, we present the numerical results on function approximation. We will detail the approximation technique in Remark 5, while the details on the explicit construction of prolongation matrices for the computation of the error will be postponed to Appendix C.
Let us consider the following smooth functions in Q = (0, 1) 3 that exhibit singularities at the origin x = (0, 0, 0): where α > 0 defines the strength of the singularity and m( 3 ) is chosen to ensure zero values of the function on Γ . Note that the function m(x) does not affect the singularity at the origin and can be represented with tensor ranks bounded from above by 3 for QTT and TQTT formats and by 9 for QT3 format.
Recall that by I we denote the Lagrange interpolation operator on the uniform tensor mesh T : In practice, u qtd will be an approximation of I u obtained by applying to A I u the exponential sums representation (see Remark 5) and by interpolating on a staggered grid (see Remark 6). We introduce the rank-truncated representation of u qtd , qtd ∈ {qtt, tqtt, qt3} based on the rounding procedure: (A u qtd , δ) , where round qtd is a rounding operation that aims at reducing the numerical qtd-rank of A u qtd with the relative Euclidean error threshold δ. The rounding procedure is based on a sequence of QR and SVD decompositions, see [65,Alg. 2] for TT (covers QTT and QT3 cases) and [18,Alg. 1] for two-level QTT-Tucker (covers the TQTT case with minor modifications).
For given u ,δ qtd , we approximate the error ε in the seminorm | · | H 1 (Q) : by using the respective quantized tensor approximation of u obtained on an equispaced mesh of axiparallel cubes with L := 30 levels of binary refinement of Q = (0, 1) 3 : Here In Fig. 6, we present the convergence plots of the relative error ε defined in (34) for δ = 10 −10 versus the background mesh levels for different α and for different quantized tensor formats. In all the cases, we observe empirical convergence in close correspondence with the rate This can be anticipated from classical FE interpolation error bounds on an equispaced, cartesian mesh in Q, for functions [x → |x| α ] in three spatial dimensions. Let us first fix α = 3/2. In Fig. 7, we present ε versus the effective number of degrees of freedom N dof for three different tensor formats. On each gray dotted line we plot the error ε for one fixed and for various values of δ. The envelopes of the computed errors with respect to N dof are highlighted with large empty markers.
In Fig. 8, we depict ε versus N dof for α = 3/2, 3/4, 1/3, and 1/4 obtained as envelopes of the set of points obtained for different δ (see Fig. 7 for α = 3/2). By plotting log 10 log 2 ε −1 against log 10 N dof , we numerically estimate the constant κ in the empirical exponential rate of convergence for some positive constants b and C. Indeed, by first applying log 2 to both sides of (35), we arrive at log 2 ε −1 =bN Assuming log 2 C is small compared with N 1/κ dof and taking log 10 of both sides, we obtain log 10 log 2 ε −1 ≈ κ −1 log 10 N dof + log 10b .
In all the numerical examples considered, we observe κ < 6, i.e., higher convergence rates than those predicted by our quantized tensor rank bounds. We also observe lower rates of convergence than those of hp-FE approximations of corner singularities in three spatial dimensions (see (43)), i.e., we find κ > 4. Figures 7 and 8 illustrate the fact that in the range of considered, the transposed order QTT representation requires more degrees of freedom to achieve a given accuracy ε than the two other formats even though it has empirical convergence (35) with slightly smaller values of κ. Among the tensor formats and for the examples considered, the TQTT format requires the smallest number of degrees of freedom to achieve a prescribed accuracy ε.
Remark 5 (Approximation of singular functions by exponential sums) To numerically evaluate the relative errors ε for all functions under consideration we used the following procedure. For each background mesh level , we approximated the function using the exponential sums representation. Specifically, we obtain the quantized tensor representations by applying the trapezoidal quadrature rule on a uniform mesh Fig. 7 Each gray dotted line represents dependence of the estimated relative seminorm | · | H 1 (D) error values ε on the rounding parameter δ for fixed as suggested by (34). Empty markers represent convex envelope of the gray dotted lines. The limits for both axes coincide for each of the plots for different values of β. A quadrature rule on a uniform mesh applied to (36) leads to an approximate, separated representation: where |x| = (x 2 1 + x 2 2 + x 2 3 ) 1/2 . Note that (37) only gives us separation of physical variables, while the exponential convergence for Gaussian functions in the QTT format was shown in [19].
The size of the integration interval and the number of points was tuned separately for each beta to ensure the desired accuracy. In this way, an approximation for |x| β , with β ∈ (0, 2), is found by first approximating the radial function |x| β−2 since β − 2 < 0 and (37) is applicable, and subsequently by multiplying this function by |x| 2 = x 2 1 + x 2 2 + x 2 3 , which has bounded TT ranks: |x| β = |x| 2 |x| β−2 for β ∈ (0, 2). This allows us to avoid using cross approximation techniques which may experience stability issues at high accuracies (using exponential sums, we obtain approximations with relative accuracy 10 −11 in L 2 norm). Note that the exponential sums approach can be applied to any of the considered TT formats: QTT, QT3 and TQTT. In this section, for the QTT-, QT3-formatted arrays and for intermediate computations in TQTT we utilized the ttpy library. 2 Remark 6 (Interpolation on staggered grid) To conveniently assemble |x| β−2 for β ∈ (0, 2) using exponential sums, while avoiding evaluation at the origin where the function has a singularity, we approximate each u(x i,j,k ) as an average of the neighboring points on a staggered grid. Let h = 2 − and denote by the nodes on the staggered grid. Then, for each x i,j,k , the set of neighboring points to x i,j,k on the staggered grid is N i,j,k = {x ∈ N : |x − x i,j,k | ≤ h /2}. We then approximate where #N i,j,k is the number of points of N i,j,k -equal to 8 except for points that lie on ∂Ω \ Γ . We need therefore function evaluations only in the points of a mesh shifted by h /2 with respect to the original mesh T , and avoid evaluations at the singularity.
After u qtd is accurately approximated for every background mesh level using exponential sums, we reduce the number of parameters in the corresponding quantized tensor representation qtd of A u qtd by using round qtd .
To solve the problem, we approximate the eigenvectors corresponding to the smallest eigenvalues in the TQTT format that yields the smallest amount of degrees of freedom for a given error (compared with the QTT and QT3 formats) according to Figs. 6 and 8. Note that due to the extremely refined underlying background meshes with 2 3 internal equispaced points, the stiffness matrix D becomes severely illconditioned (its condition number scales as h −2 , i.e., it grows exponentially in ). Besides, there arises an effect of ill-conditioning for large connected purely with the structure of tensor decompositions, see [5]. Therefore, in order to overcome the effect of algebraic and representation ill-conditioning and to accurately approximate the eigenvalues and corresponding eigenvectors of (38), particular attention has to be devoted to technical details of the computation. The overall procedure-based on the preconditioned gradient descent method and on the Rayleigh-Ritz procedure-is summarized in Algorithm 1. In the algorithm, we utilize "derivative-free" formulas [69] (that avoid multiplications by D , see Algorithm 1, line 8) for calculating the N ev × N ev matrix F given by 3 Here, linear operators are mappings A : R 2 ×2 ×2 → R 2 ×2 ×2 given as 6-dimensional arrays such that the action on u ∈ R 2 ×2 ×2 is defined by A i,j,k,p,q,r u p,q,r , (i,j,k) ∈ {1, . . . , 2 } 3 .

Algorithm 1
Block eigensolver in TQTT format based on derivative-free formulas. The algorithm is formulated for three-dimensional arrays, implying that all the operations are performed within the TQTT format.
where u ,k α are three-dimensional arrays represented in the TQTT format that approximate u α on the kth step of the iterative process and ·, · denotes scalar products of three-dimensional arrays: To solve the screened Poisson's equations arising in Algorithm 1, we utilize the algorithm proposed in [68], which is based on the alternating direction implicit method, adapted to tensor formats with the help of rank truncation after each iteration. We also note that in this solver, the stiffness matrix D is never formed explicitly in the quantized format, and each ADI step is performed in a derivative-free way. This allows us to avoid the ill-conditioning of QTT [5] and approximate the solution without stability issues.
Let λ n,l,m (δ), n = 1, 2, 3 (N ev = 14) be the eigenvalues obtained by using Algorithm 1 with a tolerance parameter δ and sorted by their quantum numbers. Let us To each λ n,l , we associate a number of degrees of freedom, which is averaged in m by analogy with (40). For every , select the parameters δ as the largest numbers satisfying where we chose δ ref = 10 −10 and where the constants c n,l satisfy c n,l > 1 (the practical choice is c n,l = 1.01). In Fig. 9, we present the errors ε = |λ n,l (δ ) − λ n | |λ n | , in eigenvalues λ n , n = 1, 2, 3 with respect to the effective number of degrees of freedom for the eigenvalue problem (38). Note that in this section, the implementation is done using the open source library TT-Toolbox 4 , which contains the implementation of the two-level QTT Tucker format [18]. In three space dimensions, this format is equivalent to the TQTT format with negligible overhead. 5

Conclusion
We considered several formats of quantized tensor-train decompositions and proved tensor rank bounds for the approximation-with a prescribed error ε ∈ (0, 1) in 4 https://github.com/oseledets/TT-Toolbox 5 In the two-level QTT Tucker format, the Tucker core of size R 1 × R 2 × R 3 is additionally decomposed using the TT decomposition, which leads to TT cores of sizes R 1 × R 1 , R 1 × R 2 × R 3 , R 3 × R 3 . So, compared with TQTT, the two-level QTT Tucker format leads to the storage of O(R 2 1 + R 2 3 ) additional degrees of freedom. H 1 (Q)-of several classes of Gevrey-smooth functions in the unit cube Q = (0, 1) 3 , with one point singularity situated at the origin. In particular, we considered singularities from Gevrey-type and analytic function spaces with regularity quantified by corresponding derivative bounds in weighted Sobolev norms, with radial weight. For these singularities, we extended the hp approximation error analysis in [73,76] to Gevrey-regular solutions with an isolated point singularity.
We then addressed approximation rate bounds in three concrete quantized tensor formats: the quantized tensor train (QTT), the transposed quantized TT (QT3) and the Tucker quantized TT (TQTT) format. Our theoretical TT rank analysis indicated that the tensor ranks and number of degrees of freedom necessary to achieve a prescribed accuracy ε ∈ (0, 1) in norm H 1 (Q) in these format might depend on the format adopted in the quantized approximation (as no lower bounds were shown, these conclusions might be an artifact of our proofs). Numerical results, however, for several model singular functions confirmed the relative rank bounds for the three mentioned formats. These results point the way to QTT-structured solvers for electron structure problems and for other PDE models where solutions exhibit isolated point singularities; for example, continua with point defects, nonlinear Schrödinger and parabolic PDEs with blowup, to name but a few. Format-adaptive, quantized approximations as were recently proposed in [7,8,61] might result in further quantitative improvement of TT ranks for the presently considered examples.
While our analysis focused only on functions with singular support consisting of one isolated point, we emphasize that corresponding rank bounds are obtained for functions whose singular support consists of a finite number of (well-separated) isolated points; the present results imply the same rank bounds as shown here also for such functions, albeit with the constants in the estimates strongly depending on the separation of the singular supports. With further analysis, the present results extend to other forms of singularities, such as line and face singularities. The details on this shall be reported in [57].

Appendix A. hp approximation in weighted Gevrey classes
We prove, in this section, the exponential convergence of the hp approximations to functions in the weighted Gevrey class J γ (Q; C, A, d) for C, A > 0, γ > 3/2, d ≥ 1. Specifically, this corresponds to functions u ∈ H 1 (Q) such that We recall that the hp space is defined as The central (novel) result of this section is the existence-for Gevrey-regular functions in J γ (Q)-of an exponentially convergent, H 1 (Q) conforming hp-projector on 1-irregular geometric meshes of hexahedra, as stated in the following proposition.
provided the uniform polynomial degree is p ≥ c 0 d for a sufficiently large constant c 0 > 0 (depending on the constant A in Eq. (41) and on d) which is independent of . The constants C hp , b hp depend on the constants C, A, and d in J γ (Q). In terms The rest of the section will be devoted to an overview of the construction of the conforming projector ,p hp . This projector has already been exploited and analyzed in detail, e.g., in [73,74]; here, we wish to sketch its construction for the sake of selfcontainedness and to provide the necessary detail of the treatment of non-analytic Gevrey-type estimates (i.e., of the cases where d > 1), which requires some minor modification with respect to the setting of [73,74]. For positive integers p and s such that 1 ≤ s ≤ p, we write Ψ p,s = (p − s)!/(p + s)!.

A.1 Discontinuous projector
We start by introducing a nonconforming projector.

A.1.1 Local projector
We denote the reference interval by I = (−1, 1) and the reference cube by K = (−1, 1) 3 . We write also H 2 mix (K) = H 2 (I ) ⊗ H 2 (I ) ⊗ H 2 (I ), where ⊗ denotes the Hilbertian tensor product. Let p ≥ 3: as constructed in [17,Appendix A], there exist univariate projectors π p : H 2 (I ) → P p (I ) such that see [73,Lemma 4.1] (the projector π p is denoted π p,2 there). Then, the Hilbertian tensor product projector given by has the following property. For all K ∈ G , we introduce the affine transformation from the reference element to K K : K → K such that K (K) = K; it follows that for v defined on K such that v • K ∈ H 2 mix (K) we can define the local projector on K so that The projector K p is continuous across regular matching faces.
Lemma 7 [73,Lemma 4.2] Let K 1 , K 2 be two axiparallel cubes that share one regular face F (i.e., F is a full face of both K 1 and K 2 ). Then, for v ∈ H 6 (int(K 1 ∪ K 2 )) the piecewise polynomial is continuous across F . (44) and (45), if a function v on K such that v • K ∈ H 2 mix (K) vanishes on a face F ⊂ ∂K, then we also have

A.1.2 Globally discontinuous hp projector
Starting from the local, elementwise projector (46), a global, discontinuous projection operator ,p hp,disc is defined in the usual way: with the nonconforming hp-space Proof The proof follows along the lines of the proof of [75,Proposition 5.13]. Denote . We start by considering K ∈ L j for j ≥ 1 and write d K = dist (K, (0, 0, 0)). By Lemma 6, scaling inequalities (see [75,Equations (5.26)-(5.31)]), and the regularity of u (see (41)), Then, using the fact that for sufficiently large s and c = 2A u + 1, [20,Equation (42)], choosing s = (p/c) 1/d , with c > 1 sufficiently large, and summing over all mesh layers not touching the origin ("interior mesh layers"), we obtain that there exist C 1 , b 1 > 0 such that for every ≥ 1 holds We now consider the element K ∈ L 0 , i.e., K = (0, 2 − ) 3 . By Hardy's inequality and choosing γ > 1, Finally, the dimension of the Q p (K) space in each element K ∈ G is given by (1 + p) 3 ; since each non-terminal mesh layer L j , j > 0, contains 7 elements, we have that dim(X ,p hp,disc ) = (1 + p) 3 (1 + 7 ). The observation that p d concludes the proof.

A.2 Conforming hp approximation
A conforming hp approximation is obtained by locally lifting the polynomial face jumps of the discontinuous, piecewise polynomial approximation. Their construction is detailed in [73,Section 5.3].

A.2.1 Edge and face liftings
Since our discontinuous interpolant is the same as in [73], apart from the nonzero constant in L 0 (see (47) and [73,Equation (4.10)]), the construction of the polynomial face jump liftings can be replicated verbatim as in [73]. We recall it here briefly, referring the reader to the aforementioned [73, Section 5.3] for the details.
We start by considering the interface between two mesh levels L k and L k+1 , k ∈ N. We introduce a local coordinate systemx,ŷ,ẑ and label the faces and edges belonging to the interface as F i , i = 1, 2, 3 and E i , i = 1, . . . , 9, respectively, see Fig. 10. Furthermore, we denote by h E the maximum length of all edges E i . We refer to Fig. 10 for the precise numbering of edges and faces and for the location of the local system of coordinates. Given two neighboring elements K a and K b with on f ab is given by where n K a (resp. n K b ) is the normal pointing outwards from element K a (resp. K b ). Consider face F 1 of Fig. 10: we define the jump of the discontinuous interpolant on this face as where f 1j are the four parts of the face F 1 , see Fig. 10. The jumps on the other faces are defined similarly. The edge jump, e.g., on E 1 , is then defined as Fig. 10 Separation of mesh levels L k (elements moved to the left) and L k+1 , with the interfaces F 1 , F 2 , F 3 and edges E 1 , . . . , E 9 marked. The local system of coordinates is given byx,ŷ,ẑ is also represented (witĥ y pointing upwards from the page) . Let n denote the normal on face F 1 pointing outwards from K 010,k+1 ; the lifting of the jump on edge E 1 is given by (49) After having defined the other edge liftings L E i , i = 2, . . . , 9, in the same way, we can introduce the full edge lifting operator We now introduce the face lifting operator for the face F 1 , the other liftings being derived in the same way. We have (50) where n is again the normal on face F 1 pointing outwards from K 010,k+1 . Then, the global lifting L k between mesh levels L k and L k+1 is the sum of the local liftings on the three interfaces: Note that the lifting thus defined has support only in the elements belonging to mesh level L k . We now turn to the terminal layer, i.e., to the jumps of ,p hp,disc u between the element K 000,0 = (0, 2 − ) 3 and the elements of L 1 . The (three) faces belonging to the interface are all regular, but ,p hp,disc u is defined as a constant in K 000,0 , see (47). One has to lift the nodal jumps at all the nodes of K 000,0 except the origin. Then, the same procedure as for the other mesh layers (applied to the nodally lifted polynomial) gives a lifting operator L 0 .
The full lifting operator is thus given by the sum of the local operators, as with all L k constructed as in (51). Such a lifting permits to obtain a conforming projector into X ,p hp , with approximation error bounded by a multiple of the approximation error of the discontinuous operator ,p hp,disc , as stated in the next proposition, that is proven in [73].  H 1 (Q) and there exists C > 0 independent of p such that hp,disc u. The exponential convergence of the conforming approximation, stated in Proposition 3, is a direct consequence of the last results. (42) follows from Proposition 4 and Lemma 8, once the algebraic term in p of inequality (48) has been absorbed in the exponential by a change of constants.

A.3 Combination of patches
We conclude this section by considering the approximation in a domain which contains the singular point in its interior. Let then R = (−1, 1) 3 . The definition of the weighted space follows directly from the definition of the spaces in Q, by keeping the weight r = |x| to be the distance from the origin.
The construction of the graded mesh is done by decomposing R into eight subcubes of unitary edge and by collecting the elements of the sub-meshes (called here "patches") obtained by symmetry from G . The projector ,p hp in R can also be straightforwardly constructed by combining local projectors obtained by symmetry; we show that it is continuous on interpatch faces, hence conforming on the whole cube R.
We detail the construction for two patches; the rest follows by iterating this argument. Specifically, we consider the two cubes and introduce the reflection operator Note that (ψ ± ) −1 = ψ ± . Then, the mesh on the domain Q ± = Q + ∪ Q − is given by see Fig. 11. The projection operator for functions v ∈ J ∞ γ (Q ± ) can be easily constructed by reflection The operator thus obtained is continuous hence conforming, as discussed in the next lemma.

Proof
,p,± hp is continuous in the sub-patches Q + and Q − . It remains to check the continuity across the interpatch interface F ± = {0} × (0, 1) 2 . By construction, all elemental faces belonging to the interface are regular, hence, by Lemma 7, the discontinuous projector ,p hp,disc is continuous across these faces. We consider the error contribution from interior mesh layers, i.e., from all elements in L j , j > 0. For all faces F in interior mesh layers which are situated Fig. 11 The mesh patches G ,− on Q − and G on Q + . An edge E belonging to the interpatch interface is highlighted perpendicular to F ± , we have We now consider any edge E belonging to F ± and separating the mesh levels L k and L k+1 , see Fig. 11 for an example. By the continuity of the discontinuous projector across regular faces where ,p,− hp,disc is the discontinuous projector in patch G ,− . Therefore, from definitions (49), (50), and (51), we conclude that the projection operator is continuous across interior mesh layers L k , k > 0.
When dealing with the terminal layer L 0 , we note that the discontinuous projector is constant hence continuous. The nodal liftings are continuous by the symmetry of their construction; the edge liftings are then continuous by the same argument as in interior mesh layers, and this gives the continuity between patches.
Finally, (53) follows from the application of the corresponding approximation results in each patch.
Note that the only modification with respect to Definition 3 is the convention r 0 = 8. -A admits a patchwise, transposed order QTT decomposition if A m,i,j,k = U 1 m,: (i 1 j 1 k 1 ) · · · U (i j k ) (56) with cores as in Definition 4 and with the restriction on the ranks r 0 = 8, r = 1. -A admits a patchwise, Tucker QTT decomposition if where, clearly, the Tucker core is now a four-dimensional array of size 8 × R 1 × R 2 × R 3 .
Let T ,R be the tensor product mesh on R given by degrees of freedom, with a positive constant C independent of ε and κ = ⎧ ⎪ ⎨ ⎪ ⎩ 4d + 3 for patchwise classic QTT , 6d + 1 for patchwise transposed order QTT , 3d + 3 for patchwise Tucker QTT .
Proof Here, we retrace the steps of the proofs of Lemmas 2, 3, and 5, generalizing them to the multipatch case.
By multiplying U 0 (m) and U 1 (i 1 j 1 k 1 ) for all m = 0, . . . , 7 and i 1 , j 1 , k 1 ∈ {0, 1}, we obtain a representation of the form (56). Patchwise Tucker QTT We Tucker-decompose the tensor V qtd , thus obtaining, by the same arguments that we used for (32), where R T ≤ C d+1 . Then, by contracting the core G and the factor Z over the index β 0 and by deriving the existence of the block QTT decomposition of the Tucker factors U , V , W as in (33), we obtain the existence of a representation of V qtd of the form (57).

Remark 9
In Proposition 5, we consider the approximation of functions in the cube R = (−1, 1) for ease of notation. Nonetheless, the argument and the result extend, without major modification, to R = (−a 1 , b 1 )×(−a 2 , b 2 )×(−a 3 , b 3 ), with a i , b i > 0, i = 1, 2, 3, and with a point singularity at the origin. This implies, by translation, that given a cube of fixed size, we can obtain bounds on patchwise quantized tensor representations that are uniform in the location of the singularity.

Appendix C. QTT representation of prolongation matrices
In order to evaluate the error ε , we need a tensor of ∇u ,δ qtd evaluated on the background mesh with L levels of refinement and have it represented using the respective tensor decomposition without accessing all its elements. This can be implemented as a multiplication by the prolongation matrices in the respective tensor format. To introduce the prolongation matrices, we start by considering the one-dimensional piecewise linear space on the background mesh with levels (recall that I j = (2 − j, 2 − (j + 1))) X 1d,0 = {v ∈ H 1 ((0, 1)) : v(1) = 0 and v | I j ∈ P 1 (I j ), j = 0, . . . , 2 − 1}.
We now show, in Proposition 6 below, that P ( →L) p.l.