Rank structured approximation method for quasi--periodic elliptic problems

We consider an iteration method for solving an elliptic type boundary value problem $\mathcal{A} u=f$, where a positive definite operator $\mathcal{A}$ is generated by a quasi--periodic structure with rapidly changing coefficients (typical period is characterized by a small parameter $\epsilon$) . The method is based on using a simpler operator $\mathcal{A}_0$ (inversion of $\mathcal{A}_0$ is much simpler than inversion of $\mathcal{A}$), which can be viewed as a preconditioner for $\mathcal{A}$. We prove contraction of the iteration method and establish explicit estimates of the contraction factor $q$. Certainly the value of $q$ depends on the difference between $\mathcal{A}$ and $\mathcal{A}_0$. For typical quasi--periodic structures, we establish simple relations that suggest an optimal $\mathcal{A}_0$ (in a selected set of"simple"structures) and compute the corresponding contraction factor. Further, this allows us to deduce fully computable two--sided a posteriori estimates able to control numerical solutions on any iteration. The method is especially efficient if the coefficients of $\mathcal{A}$ admit low rank representations and algebraic operations are performed in tensor structured formats. Under moderate assumptions the storage and solution complexity of our approach depends only weakly (merely linear-logarithmically) on the frequency parameter $1/\epsilon$, providing the FEM approximation of the order of $O(\epsilon^{1+p})$, $p>0$.


Introduction
Problems with periodic and quasi-periodic structures arise in various natural sciences models and technical applications. Quantitative analysis of such problems requires special methods oriented towards their specific features. For perfectly periodic structures, efficient methods are developed within the framework of the homogenization theory (see, e.g., [1,3,8] and other literature cited therein). However, classical homogenization methods cover only one class of problems (all cells are self-similar and the amount of cells is very large). In this paper, we use a different idea and suggest another modus operandi for quantitative analysis of boundary value problems with periodic and quasi-periodic coefficients. It generates approximations converging (in the energy space) to the exact solution and provides guaranteed and computable error estimates. The approach is applicable to (see, e.g., Figures 1, 2) (i) periodic structures, in which the amount of cell is considerable (e.g., -) but not large enough to neglect the error generated by the respective homogenized model; (ii) quasi-periodic structures that contain cells with defects and deformations; (iii) multi-periodic structures where the coefficients reflect the combined effect of several functions with different periodicity. In general terms, the idea of the method is as follows. We consider the problem P: where V is a reflexive Banach space with the norm ‖ ⋅ ‖ V , V * is the space conjugate to V (the respective duality pairing is denoted by ⟨v * , v⟩), and A : V → V * is a bounded linear operator. It is assumed that the operator A is positive definite and invertible, so that problem (1.1) is well posed. However, P is viewed as a very difficult problem because A is generated by a complicated physical structure, which may contain a huge amount of details. Therefore, attempts to solve (1.1) numerically by standard methods may lead to enormous expenditures. Similar difficulties arise if we wish to verify the quality of a numerical solution.
Assume that the operator A is approximated by a simplified positive definite operator A ∘ and the inversion of A ∘ is much simpler than the inversion of A. By means of A ∘ , we construct an iteration method based on solving a "simple" problem P : A ∘ u ∘ = g. In other words, the method is based on the operation g → A − ∘ g. It also includes the operation v → Av, which can be performed very efficiently by tensor-type decomposition methods provided that physical structures generated A have low-rank representations. We prove that iterations generate a sequence of functions converging to the exact solution of (1.1) with a geometrical rate. Furthermore, we deduce explicitly computable and guaranteed a posteriori error estimates adapted to this class of problems. They evaluate the accuracy of approximations computed on each step of the iteration algorithm. These estimates also use only inversion of A ∘ and operations of the type v → Av. In the iteration methods and error estimates inversion of the operator A is avoided.
In this paper, we consider one class of problems associated with divergent type elliptic equations where A = Q * ΛQ and A ∘ = Q * Λ ∘ Q. Here Λ : Y → Y is a bounded operator induced by a complicated quasi-periodic structure while Q : V → Y and Q * : Y → V * are conjugate operators, i.e., (y, Qw) = ⟨Q * y, w⟩ for all y ∈ Y and w ∈ V, where Y is a Hilbert space with the scalar product (⋅, ⋅) and the norm ‖ ⋅ ‖. The operators Q and Q * are induced by differential operators or certain finite-dimensional approximations of them. Henceforth, it is assumed that f ∈ V, where V is a Hilbert space with the scalar product (⋅, ⋅) V . This space is intermediate between V and V * , i.e., V ∈ V ∈ V * .
Then, the structural operators Λ and Λ ∘ are spectrally equivalent: c (Λ ∘ y, y) ≤ (Λy, y) ≤ c (Λ ∘ y, y), (1.2) where the constants are the minimal and maximal eigenvalues of the generalized spectral problem Λy − μΛ ∘ y = . Obviously, they satisfy the estimates c ≥ λ ⊖ /λ ∘ ⊕ and c ≤ λ ⊕ /λ ∘ ⊖ (which may be rather coarse). Concerning the operator Q, we assume that there exists a positive constant c such that Generalized solutions of the problems P and P are defined by the variational identities In Section 2, we show that a sequence {u k } converging to u in V can be constructed by solving problems (1.4) with specially constructed right-hand sides f k generated by the residual of (1.3). In proving convergence, the key issue is analysis of the spectral radius of the operator (1.5) and selection of such relaxation parameter ρ that provides the best convergence rate. Moreover, iteration procedures of such a type become contracting if the iteration parameter is properly selected. This fact is often used in proving analytical results (see, e.g., [20], where classical results on existence and uniqueness of a variational inequality are established by contraction arguments). Also, these ideas were used in the construction of various numerical methods (see, e.g., [7]). However, achieving our goals requires more than the fact of contraction. We need explicit and realistic estimates of the contraction factor (which are used in error analysis) and a practical method of finding Λ ∘ with minimal q. The latter task leads to a special optimization problem that defines the most efficient "simplified" operator Λ ∘ among a certain class of "admissible" operators. This question is studied in Section 3. In general, Λ and Λ ∘ can be induced by scalar, vector, and tensors functions. We show that selection of the optimal structural operator Λ ∘ is reduced to a special interpolation type problem, which is purely algebraical and does not require solving a differential problem (therefore a suitable Λ ∘ can be found a priori). We discuss several examples and suggest the corresponding optimal (or quasi-optimal) Λ ∘ , which guarantees convergence of the iteration sequence with explicitly known contraction factor. Now, it is worth discussing the main differences between our approach and the classical homogenization method developed for regular periodic structures. This method operates with a homogenized boundary value problem Q * Λ H Qu H = f , where Λ H is defined by means of an auxiliary problem with periodical boundary conditions in the cell of periodicity. The respective solution u H contains an irremovable (modeling) error depending on the cell diameter ϵ. Moreover, if ϵ tends to zero, then typically u H converges to u only weakly (e.g., in L ). Getting a better convergence (e.g., in H ) requires certain corrections, which lead to other (more complicated) boundary value problems in the cell of periodicity. The respective "corrected" solution u c H also contains an error. Typically, the error is proportional to ϵ and can be neglected only if the amount of cells is very large. If our method is applied to perfectly periodical structures then setting Λ ∘ := Λ H is one possible option. In this case, the homogenized operator (defined without correction procedures) is used for a different purpose: construction of a suitable preconditioning operator. The latter operator generates numerical solutions converging to the exact solution in the energy norm (i.e., the method is free from irremovable errors) and can be applied for a rather wide range of ϵ. In addition, the theory suggests other simpler ways of selecting suitable Λ ∘ . In this context, it is interesting to know whether or not the choice Λ ∘ := Λ H always yields the minimal value of the contraction factor. In Section 3, we briefly discuss this question and present an example of that the best Λ ∘ may differ from Λ H .
In Section 4, we deduce a posteriori estimates that provide fully computable and guaranteed estimates of the distance to the exact solution u for any numerical approximation u k,h computed for an approximation subspace V h . These estimates are established by combining functional type a posteriori estimates (see [22,25,26] and references cited therein) and estimates generated by the contraction property of the iteration method (see [24,29]).
The second part of the paper is devoted to a fast solution method for the basic iteration problem (2.1). The key idea consists of using tensor-type representations for approximations, what is quite natural if both coefficients of the respective quasi-periodic structure and the right-hand side admit low-rank tensor-type representations. We notice that the amount of structures representable in terms of low-rank formats is much larger than the amount of periodic structures covered by the homogenization method. The idea of tensor-type approximations of partial differential equations traces back to [9]. In computational mechanics this method is known as the Kantorovich-Krylov (or extended Kantorovich) method. However, it is rarely used in modern numerical technologies. In part, this is due to restrictions on the shape of the domain imposed by the Kantorovich method. Henceforth, we assume that the domain Ω satisfies these restrictions, i.e., it is a tensor-type domain (e.g., rectangular) or a union of tensor-type domains. This assumption induces certain geometrical  limitations. However, they can be bypassed by such methods as coordinate transformation and domain decomposition, which are widely used in modern computational mathematics (e.g., in iso-geometric analysis).
The recent tensor numerical methods (for steady state and dynamical problems) based on the advanced nonlinear tensor approximation algorithms have been developed in the last ten years. Literature survey on the modern tensor numerical methods for multi-dimensional PDEs can be found in [13,14,16]. In the context of problems considered in the paper, we are mainly concerned with another specific feature: very complicated material structure. In this case, direct application of standard finite element methods suffers from the necessity to account huge information encompassed in coefficients (especially in multi-dimensional problems). We show that tensor-type methods allow us to reduce computations to a collection of one-dimensional problems, which can be solved very efficiently using low-rank representations with the small storage requests. Similar ideas are applied for computing a posteriori error estimates.
Section 5 discusses numerical aspects of the method and exposes several examples. Typical behavior of quasi-periodic coefficients is described by oscillation around a constant, modulated oscillation around a given smooth function, or oscillation around a piecewise constant function. Figure 1 (1D case) represents examples of highly oscillating (left) and modulated periodic coefficients (right) functions. Figure 2 (2D case) illustrates the well-separable equation coefficient obtained by a sum of step-type and uniformly oscillating functions, namely, where C = , C = , a = , and ω = .
Brought to you by | Jyväskylän Yliopisto University Authenticated Download Date | 7/18/17 1:04 PM We show that specially constructed FEM type approximations of PDEs with slightly perturbed or regularly modulated periodic coefficients on d-fold n × ⋅ ⋅ ⋅ × n tensor grids in ℝ d may lead to the discretized algebraic equations with the low Kronecker rank stiffness matrix of size n d × n d , where n = O( ϵ ) is proportional to the large frequency parameter ϵ . In this case, the rank decomposition with respect to the d spacial variables is applied, such that the discrete solution can be calculated in the low-rank separable form, which requires storage size of only O(dn) instead of O(n d ) = O( ϵ d ) complexity representations which are mandatory for the traditional FEM techniques (the latter quickly leads to the bottleneck in case of small parameter ϵ > ). The arising linear system of equations can be solved by preconditioned iteration with the simple preconditioner Λ ∘ , such that the storage and numerical costs scale almost linearly in the univariate discrete problem size n. Numerical examples in Section 5 demonstrate the stable geometric convergence of the preconditioned CG (PCG) iteration with the preconditioner Λ ∘ and confirm the low-rank approximate separable representation to the solution with respect to d spacial variables even in the case of complicated quasi-periodic coefficients.
This approach is well suited for applying the quantized-TT (QTT) tensor approximation [15] to functions discretized on large tensor grids of size proportional to the frequency parameter, i.e., n = O( ϵ ), as it was demonstrated in the previous paper [17] for the case d = . The use of tensor-structured preconditioned iteration with the adaptive QTT rank truncation may lead to the logarithmic complexity in the grid size, O(log p n), see [14,16,23] for the rank-truncated iterative methods, [4,[10][11][12] for various examples of the QTT tensor approximation to lattice structured systems, and [2] for tensor approximation of complicated functions with multiple cusps in ℝ d .
In Section 6, we conclude with the discussion on further perspectives of the presented approach for 2D and 3D elliptic PDEs with quasi-periodic coefficients.

The Iteration Method
Let v ∈ V and ρ ∈ ℝ + . Consider the problem: find u v such that Obviously, the right-hand side of (2.1) is a bounded linear functional on V, so that this problem has a unique solution u v . Thus, we have a mapping T ρ : V → V, which becomes a contraction if the parameter ρ is properly selected. Indeed, for any v and v in V, we obtain Hence Since Λ and Λ ∘ are invertible with trivial kernels, μ and y μ are an eigenvalue and the respective eigenfunction of Λy μ = μΛ ∘ y μ if and only if they are an eigenvalue and the eigenfunction of the problem ΛΛ − ∘ Λy μ = μΛy μ . This means that Hence The minimum of the expression in round brackets is attained if ρ = ρ * := c /c . For ρ = ρ * , we find that Hence T ρ is a contractive mapping with explicitly known contraction factor q * . Well-known results in the theory of fixed points (see, e.g., [29]) yield the following theorem.

Remark 2.2. From (2.3) we obtain
This relation yields a simple (but not very sharp) estimate of the contraction factor.
For further analysis, it is convenient to estimate the right-hand side of (2.3) by a different method. Let ‖ ρ ‖ ∘ denote the operator norm which shows that T ρ is a contraction provided that In applications ρ is a self-adjoint bounded operator acting in a finite-dimensional space, so that verification of this condition amounts to finding ρ which yields the respective spectral radius of ρ (see Section 4).

Selection of Λ ∘
In this section, we discuss how to select Λ ∘ in order to minimize q, which is crucial for two major aspects of quantitative analysis: convergence of the iteration method and guaranteed a posteriori estimates. We assume that V, V, and Y are spaces of functions defined in a Lipschitz bounded domain Ω (namely y(x) ∈ for a.e.
x ∈ Ω where may coincide with ℝ, ℝ d , or d×d ) and the operators Λ and Λ ∘ are generated by bounded scalar functions, matrices or tensors. In this case, where ⊙ denotes the respective product of scalar, vector, or tensor functions. In view of (2.7) and (2.8), the value of ρ should minimize the quantity sup y∈Y (Λ ∘ ρ y, ρ y)/(Λ ∘ y, y). This procedure yields the contraction factor whose computation is reduced to solving algebraic problems at a.e. x ∈ Ω, i.e., Let be a certain set of "simple" operators defined a priori (e.g., it can be a finite-dimensional set formed by piecewise constant or polynomial functions). Then, finding the best "simplified" operator amounts to solving the following problem: findΛ ∘ ∈ such that Q(Λ,Λ ∘ ) is minimal. In other words, the optimalΛ ∘ is defined by the problem Notice that (3.1) is an algebraic problem, which should be solved (analytically or numerically) before computations. The respective solutionΛ ∘ defines the best operator to be used in the iteration method (2.6) and yields the respective contraction factor. Below we discuss some particular cases, where analysis of this problem generates an optimal (or almost optimal) Λ ∘ . Problem (3.1) is explicitly solvable if Λ ∘ and Λ have a special structure, namely, where is the unit operator and a ∘ (x) and a(x) are positive bounded functions defined in Ω. Then, Minimization with respect to ρ yields the best value ρ * = h ⊖ +h ⊕ and the respective value In accordance with (3.1), the identification of the optimal simplified problem is reduced to the problem sup a ∈ J(a, a ∘ ), where is a given set of functions. We illustrate the above relations by means of several examples.  In the simplest case, we set = P , i.e., a is a constant. From (3.3) it follows that q = (a − a)/(a + a), where a := min x∈Ω a(x) and a := max x∈Ω a(x). Then ρ * = a /(a + a) and the iteration procedure (2.6) with ρ = ρ * has the form Example 3.2 (Oscillation Around a Given Function). Consider a somewhat different example. Let a(x) be a function oscillating around a certain mean function g(x) so that If g is a relatively simple function, then it is natural to set a ∘ (x) = g(x). By (3.2), we find that h ⊕ = + ϵ, h ⊖ = − ϵ, and q = ϵ. Hence the method is very efficient for small ϵ (i.e., if a oscillates around g with a relatively small amplitude). Figures 1 and 2 illustrate three examples of quasi-periodic coefficients a and respective a ∘ corresponding to the case of oscillation around a constant with smooth modulation, oscillation around a given smooth function, or oscillation around a piecewise constant function.
Since the constants c i are defined up to a common multiplier, we can without a loss of generality assume that where λ i > and satisfy (3.4). If N = , then problem (3.5) has a simple solution, which shows that the ratio It is interesting to compare these results with those generated by homogenized models in the case of perfectly periodic structures. For this purpose, we consider a simple one-dimensional problem It is easy to see that H may not generate the best piecewise constant a ∘ , which produces the smallest contraction factor in the iteration procedure (2.6).

General Estimate
Since T ρ is a contractive mapping, we can use the Ostrowski estimates (see [24,26,29]) of the distance between v ∈ V and u (the fixed point). The estimates state that This estimate cannot be directly applied because v ρ := T ρ v is generally unknown (it is the exact solution of a boundary value problem). Instead, we must use a numerical approximation v ρ (in our analysis, we impose no restrictions on the method by which the function v ρ ∈ V was constructed). Thus, the difference η ρ := v − v ρ is a known function and the quantity δ ρ = ‖η ρ ‖ ∘ is directly computable. It is easy to see that To deduce a fully computable majorant of the norm ‖ v ρ − v ρ ‖ ∘ we use the method suggested in [25,26]. First, we rewrite (2.1) in the form For any y ∈ Y and w ∈ V , we have We estimate the first term on the right-hand side of (4.4) as follows: where ‖w * ‖ = sup w∈V ⟨w * , w⟩/‖w‖ ∘ is the dual norm. Hence ). In view of (4.3), ⟨Q * y − ρf, w⟩ = , and the majorant is equal to ‖v ρ − v ρ ‖ ∘ . Thus, estimate (4.5) has no irremovable gap and a properly selected y yields a sharp upper bound of the error.

Remark 4.1.
It is not difficult to show that the last term of M ⊕ (η ρ , τ y ) can be estimated via an explicitly computable quantity provided that y has the same regularity as the true flux (see [26]). However, in our subsequent analysis Q * y − ρf = and these advanced forms of the majorant are not required. In this case, the majorant has a simpler form: It is important that the computation of the majorant M ⊕ does not require inversion of the operator Λ associated with a complicated quasi-periodic problem.
where M ⊕ is defined by (4.5), τ y := y − ρΛQv, and y is a function in Y.

Remark 4.3.
Here η ρ and δ ρ are directly computable and q(ρ) is defined in accordance with relations presented in the previous section. Hence the cost of (4.6) is mainly related to the quantity M ⊕ (η ρ , τ y ), which is an a posteriori error majorant of the functional type. The derivation of such estimates is performed by purely functional methods and does not exploit special features of approximations (e.g., Galerkin orthogonality), numerical method, and exact solution (e.g., extra regularity). Properties of the majorants are well studied (see [25,26] and the literature cited therein). Numerous tests performed for different boundary value problems have confirmed high practical efficiency of error majorants of the functional type. It was shown that M ⊕ is a guaranteed and efficient majorant of the global error and generates good indicators of local errors if y is replaced by a certain numerical reconstruction of the exact dual solution. There are many different ways to obtain suitable reconstructions with minimal expenditures (concerning this point we refer to [21] where the reader will find a systematic discussion of computational aspects in the context of various boundary value problems). Error majorants of this type can be also used for the evaluation of modeling errors (see [27,28]). Usually, the cost of a good estimate (with the efficiency index between 1 and 2) is comparable with the cost of a numerical solution. However, the proportion essentially depends on the numerical method used. For the classical FEM schemes the expenditures are maximal (because this method generates rather coarse approximations of fluxes). For the dual mixed method, finite volume method, isogeometric approximations, and other methods producing locally equilibrated fluxes, the expenditures may be two to three times smaller than for the numerical solution.

Examples
Now we shortly discuss applications of Theorem 4.2 to problems where Q and Q * are defined by the operators ∇ and div, respectively, Λ ∘ = a ∘ (x) , Λ = a(x) , x ∈ Ω, and V = ∘ H (Ω).

d = .
Let Ω = ( , ). Equation (1.1) has the form (a(x)u ὔ ) ὔ − f = . In this case, Qw = w ὔ , Q * y = −y ὔ , and (4.3) is reduced to (4.7) In order to apply Theorem 4.2, we set y = ρ(g(x) + μ), where g(x) = − ∫ x fdx and μ is a constant. Then −y ὔ − ρf = and τ = ρ(g(x) + μ) − ρav ὔ = ρ(μ + g − av ὔ ). The best constant μ is defined by minimization of M ⊕ (η ρ , τ), which has the form Since ∫ η ὔ ρ dx = , the problem is reduced to minimization of the second term and the best μ satisfies the equation Hence , and (4.6) yields the estimate where Here v and v ρ are two consequent numerical approximations (e.g., finite element approximations v k,h and v k+ ,h computed on a mesh I h ). Then Since a ∘ is a "simple" function, the integrals are easy to compute. Other integrals Brought to you by | Jyväskylän Yliopisto University Authenticated Download Date | 7/18/17 1:04 PM contain a highly oscillating coefficient a multiplied by piecewise polynomial mesh functions. If a and f have low QTT rank tensor representations [15], then the integrals can be efficiently computed by tensor-type methods discussed in [17] (see also Section 5 below). We have Notice that the quantity ε k,h is the value of the majorant, where the flux has been selected in the best way. It is not difficult to show that a ∘ (v ὔ ρ − v ὔ ) = ρ(g +μ ) − av ὔ and, therefore, In other words, this term coincides with the error of the Galerkin solution related to the simplified boundary value problem (4.7), where v = v k,h . In accordance with Section 3, we set ρ = /(h ⊖ + h ⊕ ) and find that . Now (4.8) yields easily computable lower and upper bounds of the error encompassed in v k,h :

Remark 4.4.
It is worth adding comments on convergence properties of the quantities δ k,h and ε k,h entering (4.9). If the mesh T h is fixed, then v k,h tends to the Galerkin approximation u h of problem (1.1) on this mesh. This fact follows from Theorem 2.1 applied to the case where the iterations are performed on the respective finite-dimensional space V h (considered as the space V). Then, ‖v k,h − u h ‖ ∘ ≤ q k ‖v ,h − u h ‖ and for any h the term δ k,h tends to zero with the geometric rate. The quantity ε k,h is equal to the error of the Galerkin solution to the simplified problem. It has different asymptotic properties. It mainly depends on T h , and for a given mesh it does not tend to zero when k → +∞. However, for any given v (which in our example is defined by v k,h ) this term goes to zero if h → provided that the mesh satisfies the standard regularity conditions. The problem with a ∘ is assumed to be much more regular than the problem with a. Therefore, in terms of h the approximation error ε k,h (associated with a ∘ ) will decrease faster than the analogous error in the original problem (e.g., for a ∘ = const, the term ε k,h is proportional to h).
Since both quantities δ k,h and ε k,h are explicitly known, estimate (4.9) (and other analogous estimates) contains a very useful information unavailable in the context of purely asymptotic error analysis. Using this information, we can organize the computational process in the best possible way by comparing iteration and discretization errors. In this process, rapidly converging iterations with respect to k should be continued until δ k,h > ε k,h . If δ k,h ≈ ε k,h , then further iterations on the mesh T h are unable to essentially improve the numerical solution. Instead, we should refine T h , project u k,h on the refined mesh, and use it as a starting point for a new series of iterations generated by problem (4.7). For each step, we compute the right-hand side of (4.9) and stop the process when it becomes smaller than the desired tolerance.
The computation of M ⊕ for 2D problems can be also reduced to the computation of one-dimensional integrals. Certainly, on the multidimensional case the amount of integrals is much larger. However, the basic tensor decomposition methods remain the same. Below we briefly discuss them with the paradigm of a simple case where Assume that approximations are represented in the form of series formed by one-dimensional functions ϕ Here Υ is a given function, which can be defined in different ways. In particular, we set The functions Υ kl must satisfy the usual linear independence conditions in order to guarantee unique solvability of the respective approximation problem. For any smooth function w vanishing on ∂Ω, we have Thus, ‖Q * y − ρf‖ = and we can use the simplified form of M ⊕ . In the simplest case Λ ∘ = a ∘ , where a ∘ is a constant. The best y minimizes the quantity which shows that y must satisfy the relation y = ρa∇v + a ∘ ∇η ρ . We select σ kl that defines the Galerkin approximation of this function, and we arrive at the system m k= m l= σ kl Introduce the following matrices: Notice that all coefficients are presented by one-dimensional integrals, which can be efficiently computed with the help of special (tensor-type) methods (see, e.g., [15][16][17][18][19]). It is not difficult to see that where Y = {Y klst } is the fourth-order tensor. Hence the left-hand side of the system (4.12) has the form Yσ + g ( ) ⊗ g ( ) . In the right-hand side we have the term Ω a ∘ ς ij ∂ϕ and the value of M ⊕ is obtained by (4.6), (4.10), and (4.11).

Low-Rank Solution of the Discrete Equation
We consider the following elliptic diffusion equation with quasi-periodic coefficient a(x) > (whose oscillations are characterized by the parameter ϵ): where the function f corresponds to the modified right-hand side in problem (4.3). In this case Γ = ∂Ω, Λ = a , Qw = ∇w, and Q * y = − div y.
In what follows we assume that f and a admit low-rank representation i.e., where the parameters R f and R a are called the separation rank. Then one may assume that the exact FEM solution can be well approximated by u K (x) = ∑ K j= u j (x )u j (x ), where K depends on the separation rank of f and a. In some cases this important property can be rigorously proven (e.g., for the Laplacian and other closely related operators). Similar low-rank approximations can be observed for the QTT tensor approximations (see [17]). Existence of a low-rank solution means that for some moderate K we have u K ≈ u up to the rank truncation threshold.  It is approximated by the following Galerkin problem for the low-rank representation u K : where V K is a subset of V formed by functions of the type Therefore, in terms of the general scheme exposed in the introduction, the problem P is now problem (5.2) and we solve it by iterations with the help of the simplified (preconditioned) problem where f k− depends on u K k− and a ∘ is a simple function (i.e., it is representable by a sum of terms a ∘ (x ) . . . a ∘ (x ) with very simple multipliers).
Given the right-hand side, problem (5.3) is much simpler than the initial problem and the stiffness matrix associated with (5.3) is computed much easier and has a simpler (low Kronecker rank) form that allows a rank-structured representation of its inverse.

Kronecker Product Representation of the Stiffness Matrix
where the generating univariate function a (x ) has the shape of six uniformly distributed bumps of height as shown in Figure 3, right. Figure 3, left, presents the oscillating part of a 2D coefficients function, which is a (x )a (x ). Here, the coefficient bumps are displaced on the coarse grid of size L × L in such a way that bumps occupy the × central box in each of the × cells, which compose the whole L × L lattice-type decomposition of Ω (we have L = in Figure 3, i.e., the size of the coarse grid is × , while L = in Figures 4 and 5, i.e., the size of the coarse grid is × ). Hence the axis scale , , , denotes the coarse grid in both x and x that describes the construction of coefficient in detail.
The examples of other possible shapes of the equation coefficient corresponding to the cases (i), (ii) and (iii) specified in Section 1 are presented in Figures 1 and 2.
We apply the FEM Galerkin discretization of equation (5.1) by means of tensor-product piecewise affine basis functions (instead of "linear finite elements") where φ i k are 1D finite element basis functions (say, piecewise linear hat functions).
We associate the univariate basis functions with the uniform grid {ν j }, j = , . . . , n ℓ , on [ , ] with the mesh size h = /(n ℓ + ). In this construction we have N = n n . . . n d basis functions φ i . Notice that the univariate grid size n ℓ is of the order of n ℓ = O( ϵ ) designating the total problem size N = O( ϵ d ). For ease of exposition we first consider the case d = , and further assume that the scalar diffusion coefficient a(x , x ) can be represented in the form with a small rank parameter R.
The N × N stiffness matrix is constructed by the standard mapping of the multi-index i into the N-long univariate index i representing all degrees of freedom. For instance, we use the so-called big-endian convention for d = and d = : respectively. Hence all matrices and vectors are defined on the long index i as usual, however, the special Kronecker structure allows the low-storage and low-complexity matrix vector multiplications when appropriate, i.e., when a vector also admits the low-rank Kronecker form representation. In particular, the basis function φ i is designated via the long index, i.e., φ i = φ i .
First, we consider the simplest case R = and let d = . We construct the Galerkin stiffness matrix A = [a ij ] ∈ ℝ N×N in the form of a sum of Kronecker products of small "univariate" matrices. Recall that given p × q matrix A and p × q matrix B, their Kronecker product is defined as a p p × q q matrix C via the block representation C = A ⊗ B = [a ij B], i = , . . . , p , j = , . . . , q .
We say that the Kronecker rank of the matrix A in the representation above equals . Now the elements of the Galerkin stiffness matrix take the form which leads to the rank-2 Kronecker product representation where ⊗ denotes the conventional Kronecker product of matrices. Here A = [a i j ] ∈ ℝ n ×n and A = [a i j ] ∈ ℝ n ×n denote the univariate stiffness matrices and M = [m i j ] ∈ ℝ n ×n and M = [m i j ] ∈ ℝ n ×n define the corresponding weighted mass matrices, e.g., By simple algebraic transformations (e.g., by lamping of the tri-diagonal mass matrices, which does not effect the approximation order of the FEM discretization) the matrix A can be simplified to the form where D , D are the diagonal matrices. The matrix A corresponds to the FEM discretization of the initial elliptic PDE with complicated highly oscillating coefficients. The simple choice of the spectrally equivalent preconditioner A ∘ corresponds to the operator Laplacian. In this case the representation in (5.4) is simplified to the discrete Laplacian matrix in the form of rank-Kronecker sum where I and I denote the identity matrices of the corresponding size. Here the simple tree-diagonal matrices A and A represent the FEM/FDM Laplacian in 1D. This matrix will be used in what follows as a prototype preconditioner for solving the linear system of equations The matrix A is constructed in general for the R-term separable coefficient a(x , x ) with R ≥ which leads to the rank-R Kronecker sum representation  Further enhancement of the tensor approximation can be based on the application of the quantized-TT (QTT) tensor approximation which has been already applied in [17] to the 1D equations with quasi-periodic coefficients. The power of the QTT approximation method is due to the perfect low-rank decompositions applied to the wide class of function-related tensors [15]. See [17] for a more detailed discussion and a number of numerical examples.
One can apply QTT approximations to problems with quasi-periodic coefficients, which can be described by oscillation with smooth modulation around a constant value, oscillation around a given smooth function, or oscillation around a piecewise constant function, see Figure 1 and examples in [17].
Let   (E) Furthermore, the following result holds [11]: the QTT rank for the periodic amplification of a reference function on a unit cell to a rectangular lattice is of the same order as that for the reference function. The rank of the QTT tensor representation to the 1D Galerkin FEM matrix in the case of oscillating coefficients was discussed in [4,17]. Figure 6 represents the right-hand side f (x , x ) and the respective solution for the discretization to equation (5.1) (with the coefficient depicted in Figure 3) on a × -grid, where

Numerical Test on the Rank Decomposition of u
The PCG solver for the system of equations (5.6) with the discrete Laplacian inverse as the preconditioner demonstrates robust convergence with the rate q ≪ . The next example demonstrates the rank behavior in the singular value decomposition (SVD) of a matrix representing the solution vector u ∈ ℝ n ×n to equation (5.6) with the × periodic coefficient shown in Figure 4, left. Figure 7 represents the rank behavior in the SVD decomposition of the solution in the case of the × periodic coefficient.
Comparing Figures 4 and 7 indicates that the exponential decay of the approximation error in the rank parameter is stable with respect to the size of the L × L lattice structure of the coefficient, i.e., the behavior of the singular values remains almost the same for different parameters ϵ = L .
Our iterative scheme includes only the matrix-vector multiplication with the stiffness matrix A that has the small Kronecker rank R, and the action of the preconditioner defined by the approximate inverse to the Laplacian type matrix. The latter has low Kronecker rank of order R B = O(| log ε| ) as shown above.
Given rank-vector u = u ⊗ u , the standard property of the Kronecker product matrices indicates that the matrix-vector multiplication enlarges the initial rank by the factor of 2, and similar with the action of the preconditioner. Hence each iterative step should be supplemented with certain rank truncation procedure which can be implemented adaptively to the chosen approximation threshold or fixed bound on the rank parameter.

Conclusions
We present a preconditioned iteration method for solving an elliptic type boundary value problem in ℝ d with the operator generated by a quasi-periodic structure with rapidly changing coefficients characterized by a small length parameter ϵ. We use tensor product FEM discretization that allows to approximate the stiffness matrix A in the form of a low-rank Kronecker sum. The preconditioner A ∘ is constructed based on certain averaging (homogenization) procedure of the initial equation coefficients such that the inversion of A ∘ is much simpler than the inversion of A. We prove contraction of the iteration method and establish explicit estimates of the contraction factor q < . For typical quasi-periodic structures we deduce fully computable two-sided a posteriori estimates which are able to control numerical solutions on any iteration. We apply the tensor-structured approximation which is especially efficient if the equation coefficients admit low-rank representations and algebraic operations are performed in tensor structured formats. Under moderate assumptions the storage and solution complexity of our approach depends only weakly (merely linear-logarithmically) on the frequency parameter ϵ . Numerical tests demonstrate that the FEM solution allows the accurate low-rank separable approximation which is the basic prerequisite for application of the tensor numerical methods to the problems of geometric homogenization.
The approach allows further enhancement based on the quantized-TT (QTT) tensor approximation which is the topic for future research work. Another direction is related to fully tensor structured implementation of the computable two-sided a posteriori error estimates. The interesting question arises how far the presented approach can be extended to the numerical analysis of elliptic equations with rather unstructured jumping coefficients arising in stochastic homogenization, see, e.g., [6].