Low-rank tensor approximation of singularly perturbed boundary value problems in one dimension

We derive rank bounds on the quantized tensor train (QTT) compressed approximation of singularly perturbed reaction diffusion boundary value problems in one dimension. Specifically, we show that, independently of the scale of the singular perturbation parameter, a numerical solution with accuracy 0<ε<1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\varepsilon <1$$\end{document} can be represented in the QTT format with a number of parameters that depends only polylogarithmically on ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}. In other words, QTT-compressed solutions converge exponentially fast to the exact solution, with respect to a root of the number of parameters. We also verify the rank bound estimates numerically and overcome known stability issues of the QTT-based solution of partial differential equations (PDEs) by adapting a preconditioning strategy to obtain stable schemes at all scales. We find, therefore, that the QTT-based strategy is a rapidly converging algorithm for the solution of singularly perturbed PDEs, which does not require prior knowledge on the scale of the singular perturbation and on the shape of the boundary layers.


Introduction
The solution of singularly perturbed elliptic differential equations constitutes a challenge for numerical approximation. The solutions to such problems exhibit boundary layers, whose correct resolution is crucial to the accurate approximation of the problem. Since those layers can get arbitrarily small, and the variations in the gradient of the solution can consequentially get highly concentrated in space, the accurate solution of singularly perturbed problems by, e.g., low-order finite element (FE) methods requires a computationally demanding number of degrees of freedom. For this reason, more effective methods have been introduced, such as hp-FE methods [32], see also [31,Chapter 3] and [21], where some a priori knowledge on the solution is exploited to construct numerical methods requiring a smaller computational effort. In some instances, especially in high dimension, the implementation of such methods can still be cumbersome.
In this paper, we discuss and analyze the numerical solution of one-dimensional, singularly perturbed elliptic equations in tensor-compressed format. Specifically, we formally approximate the problem using low-order, piecewise linear finite elements and compress the resulting algebraic problem (neither the right-hand side nor the matrix of this system are formed explicitly in computations) using the quantized tensor train (QTT) method [16,23]. By doing this, we obtain an approximation accuracy comparable with that of a low-order finite element on a very fine grid (the so-called background grid), but requiring only a significantly smaller number of degrees of freedom to compute and represent the solution. We also remark that the QTT-based approach to the solution of singularly perturbed problems does not require a priori knowledge of the scale of the singular perturbation, nor of the explicit form of the layers, as for, e.g., enriched spectral methods [6].

Contributions of this paper
First, we theoretically prove and numerically verify that, for all 0 < < 1 , we can obtain a QTT approximation of the solution with accuracy > 0 , and that we can represent it with O(|log | ) parameters, with = 3 . For a given accuracy, the number of parameters of the QTT approximation is independent of the singular perturbation parameter. This is the main theoretical contribution of this paper, and it is stated in Theorem 1. We remark that this is a theoretical upper bound, in numerical experiments we find that can be smaller in practice.
The second contribution of this paper is the adaptation of the preconditioner introduced in [1] to this singularly perturbed case. The straightforward application of classic solvers (DMRG [33], AMEn [8], etc.) to the QTT format tends indeed to have stability issues, which greatly limit the background grid sizes that can be used in practice. In [1], a BPX preconditioner was developed to overcome this issue. We adapt this BPX preconditioner to our case in order to obtain stable solutions for all values of the perturbation parameter 0 < < 1 . With this at hand, we are able to reach a background grid size of around 2 −50 and to accurately solve problems with = 10 −16 . We remark that such an approximation represented as a full piecewise linear FE function would require approximately 10 15 degrees of freedom, while it is easily tractable in the tensor-compressed QTT format.

Tensor compressed solution of PDEs
Historically, the first appearance of tensor decompositions dates back to F. Hitchcook in [12]. In more recent years, a wide range of tensor decompositions have appeared and have been applied to many fields of science and engineering, see [17,18]. The tensor train (TT) decomposition, specifically, was introduced in [24] as an easy to construct low-rank decomposition of high-dimensional matrices and has its roots in matrix product state representations in physics [30]. Shortly after, it was realized that low-rank tensor representations can handle certain low-dimensional partial differential equation (PDEs) that exhibit rough behavior and are computationally challenging for conventional methods, through so-called quantization. This refers to the process of reshaping low-dimensional tensors with high mode sizes into high-dimensional tensors with low mode sizes and then applying tensor decomposition. By combining quantization and the TT representation, one obtains the QTT representation.
The QTT-formatted solution of PDEs proves useful only if the tensors involved have low QTT-ranks. To theoretically analyze the rank behavior of the problem under consideration, we approximate the solution with high-order piecewise polynomials and L 2 -project the resulting approximation onto the low-order piecewise linear finite element space. We then show that the constructed FE function admits an exact QTT decomposition with low QTT-ranks (that is, QTT-ranks that grow only linearly with respect to the polynomial degree of the high-degree approximation). The strategy used here is partially different from the approach in [15,20], where the highorder piecewise polynomial was interpolated. L 2 -projections have the advantage over interpolation operators of being stable with respect to the -dependent norm we compute the error in. It is worth noting that the strategy to obtain rank bounds can be extended to the multi-dimensional case, although the analysis can be significantly more complex than in the one-dimensional case.

Structure of the paper
In Sect. 2 we introduce the singularly perturbed problem and the functional setting of the paper. Section 3 contains the main theoretical findings of this paper, i.e., the rank bounds and the error analysis of the QTT-compressed approximate solution of the singularly perturbed problem. Specifically, we show in Theorem 1 that the number of parameters of the QTT representation of the approximate solution grows only polylogarithmically with respect to the approximation error. In Sect. 4 we discuss the numerical stability of the QTT-formatted problem and propose a preconditioning strategy whose implementation details are deferred to Appendix B. Finally, in Sect. 5 we present numerical experiments to verify the theoretical results obtained in 2 Page 4 of 32 Sect. 3 and the role of preconditioning. We conclude and discuss extensions of the present work in Sect. 6.

Statement of the problem
We consider the following problem on the interval I = (0, 1) where 0 < < 1 , 1 , 0 ∈ ℝ , and where For the weak formulation of problem (1), we introduce the Sobolev spaces and The weak formulation of (1) reads then: find u ∈ H 1 D (I) such that By Lax-Milgram's theorem, for all > 0 , problem (3) is well-defined and has a unique solution. We introduce, on H 1 (I) , the -dependent energy norm We do our analysis in the energy norm just introduced. In the literature, the stronger, balanced norm is sometimes instead considered, see [22]. An analysis in this norm is out of the scope of the present paper.

Notation
We write ℕ = {1, 2, … } for the set of positive natural numbers and ℕ 0 = ℕ ∪ {0} . For p ∈ ℕ 0 and ⊂ ℝ , let ℙ p ( ) denote the space of polynomials of degree at most p defined on . We use the convention that capitalized letters are used for multidimensional arrays and non-capitalized letters are used for vectors.
Throughout, if not stated otherwise, we use the convention that C > 0 denotes a generic constant independent of the perturbation parameters and independent of the discretization. C may change value without notice.
With the word tensor we indicate multidimensional arrays: a general d-dimensional tensor is an element A ∈ ℝ n 1 ×…×n d , with n i ∈ ℕ for i = 1, … , d , and requires storage of order O(n 1 ⋅ … ⋅ n d ).

Low-rank QTT approximation
In this section, we develop the error analysis and the rank bounds for the low-rank QTT-formatted approximation of the solutions to (1). We construct the low-rank representation of solutions u to (1) by constructing a piecewise, high-order polynomial approximation to u , then reapproximating it in a low-order finite element space, and finally QTT compressing the resulting vector of coefficients.
We start by introducing the TT and QTT representations in Sect. 3.1. In Sect. 3.2 we introduce the high-order and the low-order finite element spaces. Finally, we derive rank bounds and estimate the approximation error in Sects. 3.3 and 3.4, respectively.

Quantized tensor train
We now formalize the concepts of tensor trains [24] and quantized tensor trains.
The tensors V j , seen as elements of ℝ r j−1 ×n×r j , are the TT-cores of the decomposition and {r j } d−1 j=1 are the TT-ranks of the decomposition.
We assume that r 1 , … , r d−1 = r ∈ ℕ , for ease of presentation. The required storage of the representation in (5) is of order O(nr 2 d) , as each V j can be regarded as a 3-tensor in ℝ r×n×r and there are d such elements. The storage requirement of a d-tensor A ∈ ℝ n×…×n in the TT-format is polynomial in the dimension d, provided that the TT-ranks satisfy r ≤ Cd k for some C, k > 0 independent of d, instead of the exponential dependence O(n d ) on d of the storage of a tensor represented as a d-dimensional array.
We now introduce the quantized tensor train (QTT) format. Let u ∈ ℝ 2 L , for L ∈ ℕ , and assume it is indexed by i = 0, … , 2 L − 1 . Any index i ∈ {0, … , 2 L − 1} admits a binary representation; that is, there exist i j ∈ {0, 1} , for j = 0, … , L − 1 , such that We define, using the representation in Eq. (6), the L-tensor U ∈ ℝ 2×…×2 such that Definition 2 A vector u ∈ ℝ 2 L is said to admit a QTT decomposition [16,23] if the corresponding L-tensor U ∈ ℝ 2×…×2 defined in Eq. The setting that we are interested in is when u ∈ ℝ 2 L is the coefficient vector of a finite element (FE) function on a uniform grid. The notion of QTT decompositions can be extended to include functions f ∶ ℝ → ℝ:

Example 1
The exponential function e − x , for ∈ ℝ , admits a QTT decomposition with respect to T int L with QTT-ranks r = 1 , where L ∈ ℕ . Any x ∈ T int L , can be 1} , in its binary representation. Then we obtain which is a rank-1 QTT decomposition.

Example 2
Any polynomial function M of degree p admits a QTT-representation with respect to T int L with QTT-ranks r = p + 1 . See, e.g., [26,Theorem 6] for details on the decomposition and [16] for a proof of this result. A generalization to piecewise polynomials can be found in, e.g., [14,Lemma 3.7] and is included in the present manuscript as Lemma 3.
We refer the reader to [17] for a comprehensive presentation of QTT decomposition.

Definition 4 For any collection
⊂Ī of ordered points of Ī and for a degree p ∈ ℕ 0 , we define and

High order finite element approximation
We now introduce a high-order finite element approximation result that will be essential for our main result. For all > 0 , we denote by T hp = { 0 , 1 , 2 , 3 } the collection of points defined by The following result was proved in [21] and enables us to approximate the exact solution u with a piecewise polynomial of potentially high degree. (1). There exist C, b, 0 > 0 independent of such that for every ∈ (0, 0 ) and for all p

ℙ 1 finite element methods
We introduce here, and will use in the QTT-compressed computations, a finite element (FE) method with piecewise linear basis functions. We introduce the finite dimensional low-order FE spaces The discretized version of the weak formulation in equation (3) reads: find u L ∈ V L such that The problem in Eq. (11) can be written algebraically as where and where S L ∈ ℝ 2 L ×2 L is the stiffness matrix and R L ∈ ℝ 2 L ×2 L is the matrix such that being the Lagrange basis associated with T int L . A L is called the system matrix and f L is called the load vector of the system in (11). See, e.g., [5,10] for more details on finite element methods.
Let us now introduce the ℙ 1 -FE Galerkin projection L ∶ H 1 (I) → V L defined such that, for all v ∈ H 1 (I), We also introduce the We remark that the L 2 (I) projection is stable with respect to the H 1 (I) and L 2 (I) norms, hence, there exists a positive constant C stab such that, for all 0 < < 1 and for all L ∈ ℕ,

Rank bounds for QTT approximation
In this section we establish that there exists a constant C > 0 independent of such that, for all L ∈ ℕ , L 2 L u ,L admits a QTT representation with QTT-ranks bounded by CL. A couple of auxiliary technical results for the proof of Proposition 2 are deferred to the Appendix.
Proposition 2 There exists a constant C > 0 such that, for all 0 < < 1 , for all L ∈ ℕ , and for all p ∈ ℕ , denoting u the solution to (1), under the hypotheses (2) and with 0 = 1 = 0 , and denoting u ,p the approximation of u as given in Proposition 1, . By Lemma 6, v has QTT-ranks bounded by p + 9 and, by [25,Theorem 3.3], M −1 L has QTT-ranks bounded by 5. Then, by [24,Section 4.3], their product M −1 L v has QTTranks bounded by the product of the QTT-ranks; that is, bounded by 5(p + 9) . The ℙ 1 -FE coefficient vector of L 2 L u ,p is given by M −1 L v , and hence has QTT-ranks bounded by 5(p + 9) . ◻

A priori error analysis
We now turn our focus to estimating the error ‖ L 2 L u ,p − u ‖ . Our goal (i.e., the result in Theorem 1 below) is to prove that there exists C > 0 such that, for all L ∈ ℕ and all 0 < < 1 , there holds where h = 1∕(2 L + 1) . First, by the triangle inequality and stability (15), The second term on the right hand side of the equation above can be estimated according to Proposition 1. It remains to bound the first term on the right hand side of equation (16).

Error analysis for ℙ 1 finite elements
In this section we collect, from the literature, norm estimates on u in terms of the forcing term f and on the ℙ 1 -FE Galerkin projection L and we extend known results on L 2 -norm estimates for L applied to u to related H 1 -seminorm estimates. For ease of presentation, in the following results we assume that the boundary data satisfies 0 = 1 = 0 and we remove this constraint in our main result in Theorem 1.
We now recall -dependent upper bounds for the Sobolev norms of the solution to (1), under different regularity assumptions on the right hand side. The first result considers the weak assumption of f ∈ L 2 (I) and was proved in [29].
There exists C > 0 , independent of and f, such that The following stronger result from [29] holds when f ∈ H 1 0 (I).

of and f, such that
We recall the following result on the ℙ 1 -Galerkin projection. We now introduce interpolation spaces between L 2 (I) and H 1 0 (I).
where See, e.g., [2,5] for more details on interpolation spaces. The following result proved in [19] implies that the right hand side f of Eq. (1) is contained in H 1∕2,∞ (I).
The following proposition establishes the desired error estimates on the ℙ 1 -FE solution of problem (1) with homogeneous Dirichlet boundary conditions. The L 2 -norm estimates were proved in [29] and the H 1 -seminorm estimates are novel to our current knowledge.
In particular, if f ∈ H 1∕2,∞ (I) then Proof The L 2 -norm estimates are proved in [29, Theorem A.1] and we adapt their strategy to prove the corresponding H 1 -seminorm estimates. We begin with proving Assertion 1 and Assertion 3 and then use an interpolation technique to establish Assertion 2. Assertion 1: Suppose that f ∈ L 2 (I) . The H 1 -seminorm estimates follow from Propositions 3 and 4 The H 1 -seminorm estimates follow from Proposition 4 and Lemma 1 and for any f 1 ∈ H 1 0 (I) . Let u i , for i = 0, 1 , be the solution of Eq. (1) with 0 = 1 = 0 but with right hand side equal to f i . By linearity and by the triangle inequality we have and by taking the infimum over all such decompositions of f we obtain that Since we obtain the desired estimate for the H 1 -seminorm

Error estimate for the low-rank QTT approximation
We are now in a position to prove our main novel result, Theorem 1. In particular, Theorem 1 establishes existence of low-rank QTT-approximations converging exponentially fast to the solution u of problem (1) with respect to the number of degrees of freedom and uniformly in 0 < < 1. Proof As we do not assume homogeneous boundary conditions on u , we introduce where g(x) ∶= 0 + ( 1 − 0 )x . Then v satisfies the following modified differential equation with homogeneous boundary conditions Clearly, the regularity assumptions on f from Sect. 2 also hold for f − cg . In particular, by Lemma 2 and by the assumption (2), we have that f − cg ∈ H 1∕2,∞ (I) . Hence, the last assertion of Proposition 5 is applicable to v . Thus, Moreover, L g = g , as g ∈ V L (g is affine), and L ∶ H 1 (I) → V L is a projection. Therefore, 2 Page 14 of 32 Since v = u − g , we obtain As we have already stated in Eq. (16), Then, inequality (18) and Proposition 1 yield the desired error estimate Let now v ,L be the approximation to v provided by Proposition 1, with p = L . By Proposition 2, L 2 L v ,L has a QTT representation with respect to T int L with QTT-ranks of order O(L) . Then, and g is an affine function in I, hence it has a QTT representation with respect to T int L with rank bounded by L. The QTT-rank of the sum of two QTT-formatted vectors is bounded by the sum of their ranks [24]. Therefore, L 2 L u ,L admits a QTT decomposition with respect to T int L with ranks O(L) , and the number of parameters of the QTT decomposition of L 2 L u ,L is of order O(L 3 ) . The last assertion now follows directly, as there exist � C, � b > 0 independent of L, such that ◻

QTT-formatted computations and preconditioning
As in previous sections, we consider Eq. (1) and approximate the solution in a low-order ℙ 1 finite element space. To solve the arising linear system (12), one can assemble the system A L and approximate the right-hand side f L directly in the QTT format. In particular, the system matrix and the load vector can be represented in approximate low-rank QTT formulation, through a combination of exact representations [13] and adaptive sampling [27]. The system of linear equations can then be approximately solved using well-established optimization-based algorithms such as the alternating minimal energy solver (AMEn) [8]. Nevertheless, recent studies [1,7,28] have indicated that stability issues may occur when using this approach for large L.
To overcome this problem, a BPX-type preconditioner in the QTT format was proposed in [1]. Instead of the standard left BPX preconditioner, the authors proposed a symmetric two-sided version that ensures robustness in the QTT format. The approach is applicable for a wide range of elliptic PDEs, but does not provide -robustness, so we have modified the preconditioner for our purposes. Let us briefly introduce the original preconditioner and its modification. The implementation details are presented in Appendix B. For all ∈ ℕ , we introduce the piecewise linear basis functions ̂ ,j that satisfy Let then P ,L ∈ ℝ 2 L ×2 be the matrix associated with the canonical injection from span(̂ ,j ) into span(̂L ,j ) . We introduce the symmetric preconditioner so that the preconditioned system is � C L A L � C LūL = � C L f L . This preconditioner was first introduced in [1] and is a symmetric version of the classic BPX preconditioner [4]. The modification to C L that we propose to use for singularly perturbed problems reads instead where , is chosen to ensure independence of the condition number in terms of . Specifically, we choose , = min(2 − −1 , 1) . Note that , = (1 + 2 2 2 ) −1 is used in [3] for a one-sided preconditioner, so the square root is applied for its twosided version; when ≫ 2 − and ≪ 2 − , our choice has the same behavior as the choice , = (1 + 2 2 2 ) −1∕2 . Details of the assembly of (20) and of the preconditioned system matrix C L A L C L are deferred to Sect. B.

Numerical experiments
In this section we consider the following instance of problem (1): where 0 < < 1 is the perturbation parameter. The corresponding Dirichlet-Neumann problem as in Eq. (42) then readŝ In practice, due to its size, the vector u qtt vec, is not explicitly treatable. The computation of (24) is therefore only conducted in the QTT format.
We use a DMRG solver for all algebraic equations in our numerical simulations and we denote by tol > 0 the tolerance level for the termination of the DMRG iterations, measured by the Frobenius norm of the residual. We also use the same value tol as an accuracy parameter for rounding of QTT-formatted tensors; see [24] for more details on rounding in the TT-format. The values of tol used in the computations will be specified on a case-by-case basis.

Non-preconditioned system
We first consider the non-preconditioned, direct application of the DMRG solver. We assemble the FE matrices in the QTT format and then solve the resulting system with the DMRG solver. Figure 1 displays the error e obtained for different values of and e ∶= ‖u − u qtt ‖ . Fig. 1 Error e , obtained without preconditioner, as a function of the logarithm of the number of grid points, for ∈ {10 −1 , 10 −3 , 10 −6 , 10 −9 , 10 −12 } and with tol = 10 −10 with tol = 10 −10 . We observe instabilities for all considered values of . Choosing a different value of tol did not influence the results.
We remark that we see three different regions in the behavior of the error and that the first two regions correspond to the expected theoretical rates of convergence of the non-compressed system, obtained in Proposition 5. Namely, denoting h = 1∕(2 L + 1), while for large values of L the approximation is unstable.

Preconditioned system
The numerical experiments in the previous section illustrate the need for preconditioning the QTT-compressed version of the algebraic equation (12) obtained from the ℙ 1 -FE method. In this section we present the results from the simulations of the preconditioned and modified system, as described in Sect. 4. The simulations for varying values of and tol are shown in Fig. 2. We observe no instabilities, as were present in the non-preconditioned system. We also observe the same asymptotic rates observed for the stable region of the non-preconditioned  Fig. 2a and b shows that this depends on the choice of tol . The error plot we have considered so far have all considered the variation of the error for a fixed truncation value tol . Estimate (17) of Theorem 1, though, gives a theoretical exponential rate of convergence, with constants independent of . This can be realized, in practice, by choosing adaptively the truncation parameter tol , so that, for each and each L, the truncation error is of the same order of magnitude as the finite element error. Figure 3 displays the error e as a function of the number of degrees of freedom for solutions obtained with the adaptive choice of the accuracy parameter tol outlined. Specifically, for each and for each L, we choose the biggest tol such that the error obtained is at most 10% larger than the error e obtained with tol = 10 −12 .
The same quantity is displayed in Fig. 4, as a function of two different roots of the number of parameters of the QTT representation of the solution. It appears that the exponent 1/3 in equation (17) of Theorem 1 constitutes, in this case, an underestimation of the rate of convergence. To properly study the experimental rate of convergence, we start from the ansatz that there exist C, b, > 0 such that then, assuming that | | log 2 C | | is small in comparison to bN dof ,

Hence
We plot then, in Fig. 5, log 10 | | log 2 e | | as a function of log 10 N dof , and estimate the exponent . In accordance with our previous observation, the values obtained strongly indicate that = 1∕2 for this set of computations.

Conclusions and future work
We have proved that the solutions of singularly perturbed PDEs in one dimension admit low-rank QTT decompositions and that the energy norm of the error converges exponentially fast with respect to the third root of the number of parameters of the QTT decomposition, independently of the perturbation parameter. This is confirmed in a numerical test case, where we observe slightly better rates of convergence than prescribed by the theory, with errors independent of the perturbation parameter if an adaptive truncation strategy is chosen. Furthermore, the preconditioned system is stable at all perturbation scales and allows for the resolution of the boundary layer for very small perturbation parameters.
The natural extension of this work is the analysis of singularly perturbed problems in higher physical dimensions. This will be the subject of future investigation; we remark that our strategy to obtain rank bounds through high order approximation and L 2 projection extends rather naturally to this case.

A Auxiliary results
First, we recall the following lemma proved in [14] on exact low-rank QTT-representation of piecewise polynomials. Next, we prove three novel lemmas which are used in the proof of Proposition 2.

Lemma 4
For any P ∈ ℙ p (I) , with p ∈ ℕ and I ⊂ ℝ , there exist q 1,j ∈ ℙ j (I) and q 2,j ∈ ℙ p−j (I) for j = 0, … , p such that for every x, y ∈ I with x + y ∈ I Proof We prove by induction over p ∈ ℕ that for any polynomial P of degree p there exist polynomials q k of degree p − k , for k = 0, … , p , such that the following equality holds The base case p = 1 is trivial. Suppose that the statement holds for every polynomial of degree p − 1 . Let P(x) = ∑ p j=0 a j x j be a polynomial of degree p ∈ ℕ . We split the sum as where we used that ∑ p−1 k=0 a k x k is a polynomial of degree p − 1 and hence we can apply the induction hypothesis to obtain such q k ∈ ℙ p−1−k (I) , for k = 0, … , p − 1 . Then, by the binomial theorem, we have that x k q k (y).
, which establishes Assertion 1. We continue to prove Assertion 2 and 3. Fix i ∈ {1, … , n} . Let x ∈ (−1, 1) and let y ∈ ( i − h, i + h) . Then hx + y ∈ ( i − 2h, i + 2h) . As If we denote by then we can write Consider the case y ∈ ( i − h, i ) . We prove that c 1 ∈ ℙ p+k+1 (( i − h, i )) . For all y ∈ ( i − h, i ) , we have that i −y h > 0 and so we split the integral in equation (26) into two parts For each j = 0, … , p , the first integral is independent of y and the second integral is a polynomial of degree j + k + 1 in y: The latter follows from the fact that the integrand is a polynomial of degree j + k on 0, i −y h , as we have that , the desired property. The case y ∈ ( i , i + h) is proved analogously. ◻ Lemma 6 Let 0 < < 1 , and let u be solution to (1) under the hypotheses (2). Let p ∈ ℕ and let u ,p be the approximation of u as given in Proposition be the Lagrange basis associated with V L 0 . Then the vector v ∈ ℝ 2 L such that admits a QTT decomposition with QTT-ranks bounded by p + 9.
The entries of ṽ with ṽ i ≠ v i can be modified by addition or subtraction with rank-1 QTT-vectors in order to be equal to the corresponding element of v. As is at most 4. Thus, we have constructed v ∈ ℝ 2 L as a QTT-vector with QTT-ranks bounded by p + 5.
In this case, we can apply Lemma 5 to each of the subintervals ( 0 , 1 ), ( 1 , 2 ) and ( 2 , 3 ) (as 2 − 1 > 2h , and 1 − 0 = 3 − 2 > 2h ) to obtain a piecewise polynomial of degree p + 2 . The fact that 2 − 1 > 2h follows from By applying Lemma 5 multiple times, we obtain that is a polynomial of degree p in the subintervals and a polynomial of degree p + 2 in the subintervals Thus, by Lemma 3, the vector in ℝ 2 L with elements (x i ) has QTT-ranks bounded by p + 9 . As v i = (x i ) this concludes the proof. ◻

B Assembly of preconditioner
We discuss here the details on the construction of the preconditioner. We follow the construction in [1], while introducing a modification to obtain stability independent of the perturbation parameter.ṽ For all 0 < < 1 , and for f , c ∈ L 2 (I) , we introduce the singularly perturbed problem with homogeneous Dirichlet-Neumann boundary conditions of finding v ∈ H 1 (I) such that Let then Ã L ∈ ℝ 2 L ×2 L and f L ∈ ℝ 2 L be, respectively, the matrix and the right hand side corresponding to the ℙ 1 FE discretization of (29) on the grid with nodes {j2 −L ∶ j ∈ {1, … , 2 L }}. We introduce the following notation for the preconditioned system: where the matrix C L has been introduced in (20). In order to solve the linear system B LūL = g L in the QTT format, we need to assemble B L and g L . As shown in Sects. In what follows, we use normal parentheses ( ) to indicate matrices and vectors and we use square brackets [ ] to indicate block structures with matrices or vectors as elements. Superscript T on a matrix denotes the usual matrix transposition and superscript T on a block structure (with elements being matrices or vectors) refers to transposition of the individual block elements.

B.1 Operations on TT-cores
We introduce some additional notation in order to present the explicit construction of the preconditioner. We represent TT-cores as block matrices and introduce two operations on TT-cores denoted by • and ⋈ . In Definition 1 we referred to V 1 , … , V d ∈ ℝ r×n×r as the TT-cores of A with r ∈ ℕ being the TT-rank of A and n ∈ ℕ the mode size of A. We generalize the definition mentioned above and introduce TT-cores of ranks p × q and mode sizes m × n , for p, q, m, n ∈ ℕ.
As a TT-core U of ranks p × q and of mode sizes m × n is specified by the 2-tensors U [ , ] , = 1, … , p and = 1, … , q , for ease of exposition and for convenience we can specify U in the block structure We use the representation in Eq. (31) whenever we specify TT-cores. Then, by our convention for transposition, we have that or equivalently, Definition 7 ([1, Sect. 3.2]) Let m, n, p, q ∈ ℕ and let U be a TT-core of ranks p × q and of mode sizes m × n . We define the (i, j)-slice U {i,j} ∈ ℝ p×q of U as the matrix defined by In order to combine different TT-cores, we introduce two operations • and ⋈ on TT-cores. Definition 8 ([13, Definition 2.1]) Let p, q, r ∈ ℕ and let m 1 , m 2 , n 1 , n 2 ∈ ℕ . Consider two TT-cores U and V of ranks p × r and r × q and of mode m 1 × m 2 and n 1 × n 2 , respectively. The strong Kronecker product U ⋈ V of U and V is the TTcore of rank p × q and mode size m 1 m 2 × n 1 n 2 given, in terms of matrix multiplication of slices of sizes p × r and r × q , by for all combinations i k ∈ {m 1 , … , m k } and j k ∈ {1, … , n k } with k = 1, 2.
Definition 9 ([1, Definition 2]) Let p, p � , r, r � ∈ ℕ and let m, n, k ∈ ℕ . Consider two TT-cores A and B of ranks p × p � and r × r � and of mode size m × k and k × n , respectively. The mode core product A • B of A and B is the TT-core of rank pr × p � r � and mode size m × n given, in terms of matrix multiplication of blocks of size m × k and k × n , by for all combinations of = 1, … , p , � = 1, … , p � , = 1, … , r and � = 1, … , r � .

B.2 Construction of the preconditioner
We start by introducing some building block that will be necessary for the explicit TT-decompositions of the constituents of Eq.

B.3 Application of preconditioner
The preconditioner in the previous subsection was introduced for Dirichlet-Neumann boundary value problems with homogeneous boundary conditions. We show here how to apply it to the Dirichlet-Dirichlet case numerically approximated in Sect. 5. We suppose then, for ease of exposition, that the reaction coefficient is constant c(x) ≡ c ∈ ℝ for all x ∈ I . The first step is to consider the following problem with homogeneous boundary conditions where g(x) ∶= 1 x − 0 (x − 1) , whose solution relates to the solution of (1) as v (x) = u (x) − g(x) . We also introduce the corresponding Dirichlet-Neumann problem for which we have introduced the preconditioner C L (20).
We use the well-known Sherman-Morrison formula [11] to transfer the solution of the preconditioned discretization of (42) to the case with Dirichlet-Dirichlet boundary conditions as in (41). Theorem 2 ([9]) Let B ∈ ℝ n×n be an invertible matrix and let u, v ∈ ℝ n be such that also B + uv T is invertible. Then the inverse of B + uv T is given by The following straightforward corollary suggests how to apply, in practice, the Sherman-Morrison formula for solving linear systems.
Corollary 1 ([9, Corollary 2]) Let B ∈ ℝ n×n and let u, v, y ∈ ℝ n . Suppose that B and B + uv T are invertible matrices and suppose that x 1 ∈ ℝ n satisfies Bx 1 = y and that x 2 ∈ ℝ n satisfies Bx 2 = u . Then satisfies (B + uv T )x 3 = y.
We can express A L in terms of a rescaling of the parts of Ã L (since the Dirichlet-Dirichlet and Dirichlet-Neumann cases have different mesh sizes) and of a rank-one correction term. Therefore, the solution of the singularly perturbed problem with homogeneous Dirichlet boundary conditions can be obtained using Corollary 1. In the case where c is nonconstant, then the preconditioned solution for the Dirichlet-Dirichlet case can be obtained in a similar fashion, provided that the integrals arising in L,0 are computed using the grid T int L .
Funding Open Access funding provided by ETH Zurich.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.