Low-rank approximations for parametric non-symmetric elliptic problems

Low-Rank Approximations for Parametric Non-Symmetric Elliptic Problems Tomás Chacón Rebollo 1 *, Macarena Gómez Mármol 2 and Isabel Sánchez Muñoz 3 1 Departamento EDAN & IMUS, Universidad de Sevilla, Sevilla, Spain, 2 Departamento EDAN, Universidad de Sevilla, Sevilla, Spain, 3 Departamento Matemática Aplicada I, Universidad de Sevilla, Sevilla, Spain In this study, we obtained low-rank approximations for the solution of parametric nonsymmetric elliptic partial differential equations. We proved the existence of optimal approximation subspaces that minimize the error between the solution and an approximation on this subspace, with respect to the mean parametric quadratic norm associated with any preset norm in the space of solutions. Using a low-rank tensorized decomposition, we built an expansion of approximating solutions with summands on finite-dimensional optimal subspaces and proved the strong convergence of the truncated expansion. For rank-one approximations, similar to the PGD expansion, we proved the linear convergence of the power iteration method to compute the modes of the series for data small enough. We presented some numerical results in good agreement with this theoretical analysis. Keywords: low-rank tensor approximations, non-symmetric problems, PGD mode computation, alternate least squares, power iteration algorithm 1 INTRODUCTION This study deals with the low-rank tensorized decomposition (LRTD) of parametric non-symmetric linear elliptic problems. The basic objective in model reduction is to approximate large-scale parametric systems by low-dimension systems, which are able to accurately reproduce the behavior of the original systems. This allows tackling in affordable computing times, control, optimization, and uncertainty quantification problems related to the systems modeled, among other problems requiring multiple computations of the system response. The PGD method, introduced by Ladévèze in the framework of the LATIN method (LArge Time INcrement method [22]) and extended to multidimensional problems by Ammar et al [2], has experienced an impressive development with extensive applications in engineering problems. This method is an a priori model reduction technique that provides a separate approximation of the solution of parametric PDEs. A compilation of advances in the PGD method may be found in [11]. Among the literature studying the convergence and numerical properties of the PGD, we can highlight [1], where the convergence of the PGD for linear systems of finite dimension is proved. In [8], the convergence of the PGD algorithm applied to the Laplace problem is proven, in a tensorized framework. The study [9] proves the convergence of the PGD for an optimization problem, where the functional framework is strongly convex and has a Lipschitz gradient in bounded sets. In [17], the authors prove the convergence of an algorithm similar to a PGD for the resolution of an optimization problem for a convex functional defined on a reflective Banach space. In [16], the authors prove the convergence of the PGD for multidimensional elliptic PDEs. The convergence is achieved because of the generalization of Eckart and Young’s theorem. Edited by: Traian Iliescu, Virginia Tech, United States Reviewed by: Yusif Gasimov, Azerbaijan University, Azerbaijan Francesco Ballarin, Catholic University of the Sacred Heart, Brescia, Italy *Correspondence: Tomás Chacón Rebollo [email protected] Specialty section: This article was submitted to Statistical and Computational Physics, a section of the journal Frontiers in Physics Received: 04 February 2022 Accepted: 28 March 2022 Published: 10 May 2022 Citation: Chacón Rebollo T, Gómez Mármol M and Sánchez Muñoz I (2022) Low- Rank Approximations for Parametric Non-Symmetric Elliptic Problems. Front. Phys. 10:869681. doi: 10.3389/fphy.2022.869681 Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696811 ORIGINAL RESEARCH published: 10 May 2022 doi: 10.3389/fphy.2022.869681 The present study is motivated by [5], where the authors present and analyze a generalization of the previous study [16] when operator Adepends on a parameter. The least-squares LRTD is introduced to solve parametric symmetric elliptic equations. The modes of the expansion are characterized in terms of optimal subspaces of a finite dimension that minimize the residual in the mean quadratic norm generated by the parametric elliptic operator. As a by-product, this study proves the strong convergence in the natural mean quadratic norm of the PGD expansion. A review of low-rank approximation methods (including PGD)maybefoundinthestudies[18,19]. In particular, minimal residual formulations with a freely chosen norm for Petrov–Galerkin approximations are presented. In addition, the study [10] gives an overview of numerical methods based on the greedy iterative approach for non-symmetric linear problems. In [3,4], a numerical analysis of the computation of modes for the PGD to parametric symmetric elliptic problems is reported. The nonlinear coupled system satisfied by the PGD modes is solved by the power iteration (PI) algorithm, with normalization. This method is proved to be linearly convergent, and several numerical tests in good agreement with the theoretical expectations are presented. Actually, for symmetric problems, the PI algorithm to solve the PGD modes turns out to be the adaptation of the alternating least-squares (ALS) method thoroughly used to compute tensorized low-order approximations of high-order tensors. The ALS method was used in the late 20th century within the principal components analysis (see [6,7,20,21]) and extended in [27] to the LRTD approximation of high-order tensors. Its convergence properties were subsequently analyzed by several authors; local and global convergence proofs within several approximation frameworks are given in [14,25,26]. Several generalizations were reported; as mentioned, without intending to be exhaustive, in the studies [12–15,23]. The convergence proofs within the studies [14,25,26] cannot be applied to our context as we are dealing with least-squares with respect to the parametric norm intrinsic to the elliptic operator, even for symmetric problems. Comparatively, the difference is similar to that between proving the convergence of POD or PGD approximations to elliptic PDEs. The POD spaces are optimal with respect to a user-given fixed norm, while the PGD spaces are optimal with respect to the intrinsic norm associated to the elliptic operator (see [5]). This use of the intrinsic parametric norm does not make it necessary the previous knowledge of the tensor object to be approximated (in our study, the parametric solution of the targeted PDE), as is needed in standard ALS algorithms. In the present study, we propose an LRTD to approximate the solution of parametric non-symmetric elliptic problems based upon symmetrization of the problem (Section 2). Each mode of the series is characterized in terms of optimal subspaces of finite dimension that minimize the error between the parametric solution and its approximation on this subspace, but now with respect to a preset mean quadratic norm as the mean quadratic norm associated to the operator in the non-symmetric case is not well-defined (Section 3). We prove that the truncated LRTD expansion strongly converges to the parametric solution in the mean quadratic norm (Section 4). The minimization problems to compute the rank-one optimal modes are solved by a PI algorithm (Section 5). We prove that this method is locally linearly convergent and identifies an optimal symmetrization that provides the best convergence rates of the PI algorithm, with respect to any preset mean quadratic norm (Section 6). We finally report some numerical tests for 1D convection–diffusion problems that confirm the theoretical results on the convergence of the LRTD expansion and the convergence of the PI algorithm. Moreover, the computing times required by the optimal symmetrization are compared advantageously to those required by the PGD expansion (Section 7). 2 PARAMETRIC NON-SYMMETRIC ELLIPTIC PROBLEMS Let us consider the mathematical formulation for parametric elliptic problems introduced in [5] that we shall extend to nonsymmetric problems. Let (Γ,B,μ)be a measure space, where we assume that the measure μis σ-finite. Let Hbe a separable Hilbert space endowed with the scalar product (·,·) and associated norm ·,denotebyH′the dual space of Hand by 〈·,·〉the duality pairing between H′and H. We will consider the Lebesgue space L2 μ(Γ)and the Bochner space L2 μ(Γ;H)and its dual space L2 μ(Γ;H′),denoting〈·,·〉as the duality between L2 μ(Γ;H)and L2 μ(Γ;H′). We are interested in solving the parametric family of variational problems: Find u:Γ→Hsuch that auγ  ,v;γ  〈fγ  ,v〉,∀v∈H, μ−a.e.γ∈Γ, (1) where a(·,·;γ):H×H→Ris a parameter-dependent, possibly non-symmetric, bilinear form and f(γ)∈H′is a parameterdependent continuous linear form. It is assumed that a(·,·;γ) is uniformly continuous and uniformly coercive on Hμ-a. e. γ∈Γand there exist positive constants αand βindependent of γsuch that, aw,v;γ  ≤βwv,∀w, v ∈H, μ−a.e.γ∈Γ,(2) αw2≤aw,w;γ  ,∀w∈H, μ−a.e.γ∈Γ.(3) By the Lax–Milgram theorem, problem (1) admits a unique solution μ-a. e. γ∈Γ. To treat the measurability of uwith respect to γ, let us consider the problem: Let fbe a function that belongs to L2 μ(Γ;H′)such that f(γ)= f(γ)μ-a. e. γ∈Γ, Find u∈L2 μ(Γ;H)such that  au,v () 〈f,v〉,∀v∈L2 μ(Γ;H), (4) where  a:L2 μ(Γ;H)×L2 μ(Γ;H)→Ris defined by Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696812 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric  aw,v () Γ awγ  ,v(γ);γ  dμγ  .(5) By the Lax–Milgram theorem, owing to (2) and (3), problem (4) admits a unique solution. Problems (1) and (4) are equivalent in the sense that (see [5]) uγ uγ  ,μ−a.e.γ∈Γ. We shall consider a symmetrized reformulation of the formulation (4). Let us consider a family of inner products in H,{(·,·)γ}γ∈Γwhich generate the norms · γ uniformly equivalent to the ·norm and there exist α H >0 and β H >0 such that, αHv≤vγ≤βHvfor all γ∈Γ,(6) considering the associated scalar product in L2 μ(Γ;H)and  Γ (·,·) γ dμ(γ). As  ais continuous and coercive on L2 μ(Γ;H),there exists a unique isomorphism in L2 μ(Γ;H)that we denote A,suchthat Γ Aw () γ  ,vγ  γdμγ  aw,v () ,∀w,v∈L2 μ(Γ;H).(7) Let us define a new bilinear form  b:L2 μ(Γ;H)×L2 μ(Γ;H)→Rby  bw,v () ≔ aw,Av () .(8) It is to be noted that  bis symmetric, as it can be written as  bw,v () Γ Aw () γ  ,Av () γ  γdμγ  .(9) Thus, as a consequence of (2) and (3), the form  bdefines a scalar product in L2 μ(Γ;H)and generates a norm equivalent to the usual norm is this space. Now, problem (4) is equivalent to Find u∈L2 μ(Γ;H)such that  bu,v () 〈f,Av〉,∀v∈L2 μ(Γ;H). . (10) We shall use this formulation to build optimal approximations of problem (1) on subspaces of finitedimension, when the form  a(·,·) is not symmetric. For a given integer k≥1, we denote by S k the Grassmannian variety of Hformed by its subspaces of dimension smaller than or equal to kand consider the problem: min Z∈Sk  bu−uZ,u−uZ () ,(11) where uis the solution of problem (4) and u Z is its approximation given by the Galerkin approximation of problem (10) in L2 μ(Γ;Z): Find uZ∈L2 μ(Γ;Z)such that  buZ,v () 〈f,Av〉,∀v∈L2 μ(Γ;Z). (12) Then, for any k∈N, an optimal subspace of dimension smaller than or equal to kis the best subspace of the family S k that minimizes the error between the solution uand its Galerkin approximation u Z on this subspace, with respect to the mean quadratic norm generated by  b. Theorem 4.1 of [5] proves the following result: Theorem 2.1 For any k≥1, problem (11) admits at least one solution. 3 TARGETED-NORM OPTIMAL SUBSPACES It is assumed that we give a family of inner products {(·,·)H,γ}γ∈Γon H,whichgeneratenorms· H,γ uniformly equivalent to the reference norm in H. Eventually, we may set (·,·) H,γ =(·,·). Our purpose is to determine the inner products (·,·) γ (introduced in Section 2 and which we will call (·,·) γ,+ )in such a way that the corresponding bilinear form  bdefined by (8) actually is  bw,v () Γ wγ  ,vγ  H,γdμγ . In this way, the optimal subspaces are the solution of the problem min Z∈SkΓ uγ −uZγ  2 H,γdμγ . (13) Let us consider the operators A γ,+ :H↦Hand the bilinear forms (·,·) γ,+ on H×Hdefined by aw,A γ,+v;γ  w, v () H,γ,∀w, v ∈H, (14) w, v () γ,+A−1 γ,+w, A−1 γ,+v  H,γ,∀w, v ∈H. (15) It is to be noted that A γ,+ is an isomorphism on Hand consequently (·,·) γ,+ is an inner product on H.Dueto(2) and (3),thenorms generated by these inner products are uniformly equivalent to the reference norm in H. Moreover, by (14) and (15), Aγ,+w, v  γ,+w, A−1 γ,+v  H,γaw,v;γ  . (16) Let us consider now the inner product in L2 μ(Γ;H)given by  Γ (·,·) γ,+ and the isomorphism A+:L2 μ(Γ;H)↦L2 μ(Γ;H) defined by Γ A+w () γ  ,vγ  γ,+dμγ  aw,v () ,∀w,v∈L2 μ(Γ;H). (17) Then, it holds. Lemma 1 Let A γ :H→Hbe the continuous linear operators defined by Aγw, v  γaw,v;γ  ,∀w, v ∈H, μ−a.e.γ∈Γ. (18) Then, it holds Aw () γ Aγwγ  ,∀w∈H, μ−a.e.γ∈Γ. (19) This result follows from a standard argument using the separability of space Hthat we omit for brevity. Then, by Lemma 1, we have A+w () γ Aγ,+wγ ,∀w∈L2 μ(Γ;H),μ−a.e.γ∈Γ. (20) Let us denote by  b+the bilinear form on L2 μ(Γ;H)given by Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696813 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric  b+w,v ()  aw,A+v () ,∀w,v∈L2 μ(Γ;H). (21) Then, by (20) and (14),  b+w,v () Γ awγ ,A γ,+vγ  ;γ  dμγ  Γ wγ  ,vγ  H,γdμγ . (22) As a consequence, the optimal subspaces obtained by the leastsquares problem (11) when  b b+satisfy (13). Remark 1 When the forms a (·,·;γ)are symmetric, if we choose w, v () H,γaw,v;γ  ,∀w, v ∈H, μ−a.e.γ∈Γ,(23) then A γ,+ defined by (14) is the identity operator in H. From (21), it follows that  b+ a.We,thus,recover the intrinsic norm in the symmetric case to determine the optimal subspaces. 4 A DEFLATION ALGORITHM TO APPROXIMATE THE SOLUTION Following the PGD procedure, we approximate the solution of problem (1) by a tensorized expansion with summands of rank ≤k, obtained by deflation. For all N≥1, we approximate u≃uN≔ N i1 si,with si∈L2 μ(Γ;H),(24) computed by the following algorithm: Initialization: let be u 0 =0. Iteration: assuming u i−1 , known for any i≥1, set e i−1 ≔u−u i−1 . Then, uiui−1+si,where siei−1 () W(25) with (ei−1)Wthe approximation of e i−1 given by the problem (12) on an optimal approximation subspace Wsolution of the problem (11), W≔argmin Z∈Sk  bei−1−ei−1 () Z,ei−1−ei−1 () Z  . (26) It is to be noted that this algorithm does not need to know the solution uof problem (4) since e i−1 is defined in terms of the current residual fif−Aui−1by ei−1∈L2 μ(Γ;H)such that  aei−1,v () 〈f,v〉− aui−1,v () ≔〈fi,v〉∀v∈L2 μ(Γ;H)  The convergence of the u N to uis stated in Theorem 5.3 of [5] as the form  bis an inner product in L2 μ(Γ;H),withthe generated norm equivalent to the standard one. Theorem 4.1 The truncated series u N determined by the deflation algorithm (24)–(26) satisfying lim N→∞  bu−uN,u−uN () 0. Consequently, u N strongly converges in L2 μ(Γ;H)to the solution uof problem (4). 5 RANK-ONE APPROXIMATIONS An interesting case from the application point of view arises when we consider rank-one approximations. Indeed, when k= 1, the solution of (12) u Z can be obtained as uZ(γ)φγ w, μ−a.e.γ∈Γ,for some φ∈L2 μ(Γ),w∈H. Then, the problem (11) can be written as (see [5], Sect. 6) min v∈H, ψ∈L2 μΓ () Jv,ψ  ,(27) where Jv,ψ   bu−ψ⊗v, u−ψ⊗v  . (28) Any solution of problem (27) has to verify the following conditions: Proposition 1 If (w, φ)∈H×L2 μ(Γ)is a solution of problem (27), then it is also a solution of the following coupled nonlinear problem:  aφ⊗w, φ⊗Aγv  〈f,φ⊗Aγv〉∀v∈H, (29)  aφ⊗w, ψ⊗Aγw  〈f,ψ⊗Aγw〉∀ψ∈L2 μ(Γ).(30) We omit the proof of this result for brevity; let us just mention that conditions (29) and (30) are the first-order optimality conditions of the problem (27) that take place as the functional J:H×L2 μ(Γ)→Rand is Gateaux-differentiable. It is to be noted that the PGD method corresponds to replacing A γ by the identity operator in (29)–(30). From Theorem 2.1 and Proposition 1, there exists at least a solution to problem (27) and then to problem (29)–(30). However, as functional Jis not convex, there is no warranty that it admits a unique minimum. Then, a solution to problem (29)–(30) could not be a solution to the problem (27). Relations (29)–(30) suggest an alternative way to compute the modes in the PGD expansion to solve (1). Indeed, we propose a LRTD expansion for ugiven by u≃uN≔ N i1 φi⊗wi,(31) where the modes (wi,φi)∈H×L2 μ(Γ)are recursively computed as a solution of the problem (fi≔f−Aui−1):  aφi⊗wi,φi⊗Aγv  〈fi,φi⊗Aγv〉∀v∈H,  aφi⊗wi,ψ⊗Aγwi  〈fi,ψ⊗Aγwi〉∀ψ∈L2 μ(Γ). ⎧ ⎨ ⎩(32) If this problem is solved in such a way that its solution is also a solution of min v∈H, ψ∈L2 μΓ () Jiv, ψ  , with Jiv, ψ   bu−ui−1−ψ⊗v, u−ui−1−ψ⊗v  , (33) then expansion (31) will be optimal in the sense of the expansion (24), where each mode of the series is computed in an optimal finite-dimensional subspace that minimizes the error. Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696814 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric 6 COMPUTATION OF LOW-RANK TENSORIZED DECOMPOSITION MODES In this section, we analyze the solution of the nonlinear problem (32) by the power iteration (PI) method. As the operator A γ appears in the test functions, a specific treatment is needed, in particular, to compute targeted-norm subspaces. We also introduce a simplified PI algorithm that does not need to compute A γ . Here, we report both. We focus on solving the model problem (29)–(30),whichwe assume to admit a nonzero solution. For simplicity, the notation fwill stand either for the r. h. s. of the problem (4) or for any residual f i . It is to be observed that (30) is equivalent to Γ φγ aw,A γw;γ  ψγ dμγ  Γ 〈fγ ,A γw〉ψγ dμγ  ,∀ψ∈L2 μ(Γ), and then φγ φw, γ  ≔〈fγ  ,A γw〉 aw,A γw;γ  μ−a.e.γ∈Γ. (34) Thus, problem (29)–(30) consists in Find w∈Hsuch that Γ φw, γ  2aw,A γv;γ  dμγ  Γ φw, γ  〈fγ ,A γv〉dμγ  ,∀v∈H. ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩(35) with φ(w, ·) ∈L2 μ(Γ)defined by (34). For simplicity, we shall denote by φ(w) the function φ(w,·). Let us define the operator T:H→Hthat transforms w∈H into T(w)∈Hsolution to the problem Γ φw, γ  2aTw () ,A γv;γ  dμγ  Γ φw, γ  〈fγ  ,A γv〉dμγ  ,∀v∈H. (36) Tis well-defined from (2) and (3), and a solution of (35) is a fixed point of this operator. Moreover, as A γ is linear, φλw () λ−1φw () and then Tλw () λTw () ,∀λ∈R. Thus, if (w,φ) is a solution to (29)–(30), then (λw,λ −1 φ) is also a solution to this problem. So, we propose to find a solution to problem (35) with unit norm. For that, we apply the following PI algorithm with normalization: Initialization:Given any nonzero w 0 ∈Hsuch that φ 0 =φ(w 0 )is not zero in L2 μ(Γ). Iteration: Knowing w n ∈H, the following is computed : i)~ wn+1Tw n () ∈H ii)wn+1~ wn+1 ~ wn+1∈H iii)φn+1φwn+1  ∈L2 μ(Γ). (37) The next result states that this iterative procedure is welldefined. Lemma 2 It is assumed that for some nonzero w∈H, it holds that φ(w)∈L2 μ(Γ)is not zero.Then ~ wT(w)is not zero in H and φ(~ w)is not zero in L2 μ(Γ). Proof First, by reduction to the absurd, it is assumed that T(w) = 0. From (36), we have Γ φw, γ  〈fγ  ,A γv〉dμγ 0,∀v∈H. In particular, for v=wand taking into account (34), we deduce that Γ |〈fγ  ,A γw〉|2 aw,A γw;γ  dμγ 0. Then, 〈fγ ,A γw〉0,μ−a.e.γ∈Γ, and thus, φ(w) = 0 is in contradiction with the initial hypothesis. This proves that T(w) is not zero. Second, arguing again by reduction to the absurd, it is assumed that φ(~ w)0. From (34),wehave 〈fγ  ,A γ~ w〉0,μ−a.e.γ∈Γ. Then, setting v~ win (36) and using (19), we obtain  bφw () ⊗~ w, φw () ⊗~ w  Γ aφw, γ  ~ w, φw, γ  Aγ~ w  dμγ  0. As  bis a scalar product in L2 μ(Γ;H), this implies that φw () ⊗~ wL2 μ(Γ;H)φw () L2 μ(Γ)~ w0. We have already proven that ~ w≠0. So, φ(w)hastobeequalto zero, in contradiction with the initial hypothesis. Thus, our assumption is false and φ(~ w)is not zero. This result proves that if w 0 and φ 0 are not zero, then the algorithm (37) is well-defined. 6.1 Computation of Power Iteration Algorithm for Targeted-Norm Optimal Subspaces From a practical point of view, in general, the algorithm (37) is computationally expensive. Indeed, in practice, His a space of large finite-dimension and the integral in Γis approximated by some quadrature rule with nodes {γi}M i1. The method requires the computation of Aγivfor all the elements von a basis of H and all the γ i . It is to be noted that when targeted subspaces are searched for, in the way considered in Section 3, the expression of algorithm (37) simplifies. Indeed, as a(w,A γ,+ v;γ)=(w,v) H,γ , then φw, γ  〈fγ  ,A γ,+w〉 w2 H,γ μ−a.e.γ∈Γ. (38) In addition, the problem (36) that defines the operator T simplifies to problem: Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696815 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric Γ φw, γ  2Tw () ,v () H,γdμγ  Γ φw, γ  〈fγ  ,A γ,+v〉dμγ ,∀v∈H. (39) We shall refer to the method (38)–(39) as the TN (targetednorm) method. 6.2 A Simplified Power Iteration Algorithm An approximate, but less expensive, method is derived from the observation that A γ is an isomorphism in H. If we approximate the first equation in (32) by  aφi⊗wi,φi⊗Aγ0v  〈fi,φi⊗Aγ0v〉∀v∈H, for some γ 0 ∈Γ, then this equation is equivalent to  aφi⊗wi,φi⊗v  〈fi,φi⊗v〉∀v∈H. (40) We then consider the following adaptation of the PI method (37) to compute an approximation of the solution of the optimality conditions (32): Iteration: Known w n ∈H, the following is computed: i)~ wn+1^ Tw n () ∈H ii)wn+1~ wn+1 ~ wn+1∈H iii)φn+1φwn+1  ∈L2 μ(Γ). (41) where ^ T(w)is computed by Γ φw, γ  2a^ Tw () ,v;γ  dμγ Γ φw, γ  〈fγ  ,v〉dμγ  ,∀v∈H, , (42) where φ(w,γ)isdefined by (34). We shall refer to method (41)–(42) as the STN (simplified targeted-norm) method. The difference between the STN method and the standard PGD one is only the definition of the function φthat in this case is given by φγ ≔〈fγ ,w〉 aw,w;γ  μ−a.e.γ∈Γ. 6.3 Convergence of the Power Iteration Algorithms In this section, we analyze the convergence of the PI algorithms (37) and (41) (Theorem 6.1). We prove that the method with optimal convergence rate corresponds to the operator A γ,+ introduced in Section 3, choosing all the inner products · H,γ equal to the reference inner product ·. As in [5], we shall assume that the iterates φ n remain in the exterior of some ball of positive radius, say ε>0, of L2 μ(Γ), that is, φnL2 μΓ () ≥ε,n0,1,2,/. (43) This is a working hypothesis that makes sense whenever the mode that we intend to compute is not zero, considering that the w n is normalized. From the definition of A γ and (6), it holds ~ αv≤Aγvγ≤~ βvfor all γ∈Γ, with ~ αα βH ,~ ββ αH ,(44) where αand βare given by (2) and (3), respectively, and α H and β H are defined in (6). Let us define the function: δr, s () 2λk21+2−r 1−rk2  1+2−r 1−rλk2s  s2, where k ~ β ~ α,and λ k2for method (37) β αfor method (41) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ It holds. Theorem 6.1 It is assumed that (43) holds, and Δδr,  s () <1, for some r ∈0,1 () with  s uL2 μ(Γ;H) ε. (45) Then, there exists a unique solution with norm 1 w of problem (35); the sequence {wn}n≥1computed by either method (37) or method (41) is contained in the ball B H (w,r)if w 0 ∈B H (w,r),and w−wn+1≤Δw−wn,∀n≥0. (46) As a consequence, the sequence {wn}n≥1that is defined by either method (37) or method (41) is strongly convergent to w with linear rate and the following error estimate holds: w−wn≤Δnw−w0,∀n≥1, whenever w0∈BHw, r () . (47) Proof. Let us consider at first the method (37). Let x∈B H (w,r) such that φ(x)L2 μ(Γ)≥ε. Denote ~ xT(x)which by the definition of operator Tin (36) is the solution to the problem Γ φx, γ  2a~ x, Aγv;γ  dμγ  Γ φx, γ  〈fγ ,A γv〉dμγ  ,∀v∈H. (48) We aim to estimate w−~ x ~ x. To do that, from problems (35) and (48), we obtain Γ φw, γ  2aw−~ x, Aγv;γ  dμγ  − Γ φw, γ  2−φx, γ  2  a~ x, Aγv;γ  dμγ  +Γ φw, γ  −φx, γ  〈fγ  ,A γv〉dμγ  ,∀v∈H. (49) It holds ay,A γz;γ  Aγy, Aγz  γ≤~ β2yz,∀y, z ∈H, (50) using (44). Thus, 〈fγ  ,A γv〉auγ  ,A γv;γ  ≤~ β2uγ  v,∀v∈H. (51) Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696816 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric Moreover, ay,A γy;γ  Aγy, Aγy  γ≥~ α2y2∀y∈H. (52) Setting vw−~ xin (62) and using (50)-(52), we have ~ α2φw () 2 L2 μΓ () w−~ x2≤ ~ β2φ2w () −φ2x () L1 μΓ () ~ x+φw () −φx () L2 μ(Γ)uL2 μ(Γ;H)  w−~ x. (53) To bound the second term in the r. h s. of (53),from(34),itholds φw, γ  −φx, γ  auγ  ,A γw−x () ;γ  aw,A γw;γ  +ax,A γx−w () ;γ  +ax−w, Aγw;γ  aw,A γw;γ  ax,A γx;γ auγ  ,A γx;γ  . Then, using (50) and (52), |φw, γ  −φx, γ  |≤k21+k21 x+1  uγ  w−x ≤k21+2−r 1−rk2  uγ  w−x, (54) where k~ β ~ α. In the last estimate, we have used that as x∈B H (w,r), w−r≤xand then 1 x≤1 1−r. Therefore, φw () −φx () L2 μ(Γ)≤ϕ1r () uL2 μ(Γ;H)w−x,(55) where ϕ1r () k21+2−r 1−rk2  . It is to be noted that from (34), |φz, γ  | auγ ,A γz;γ  az,A γz;γ  ≤k2uγ   z,∀z∈H. (56) Then, from (54) and (56) |φw, γ  2−φx, γ  2|≤|φw, γ  −φx, γ  ||φw, γ  +φx, γ  | ≤k2ϕ1r () 1+1 x  uγ  2w−x. Hence, φw () 2−φx () 2L1 μ(Γ)≤ϕ2r () u2 L2 μ(Γ;H)w−x,(57) where ϕ2r () k2ϕ1r ()2−r 1−r. Combining (53) with (55) and (57), we deduce w−~ x≤ uL2 μ(Γ;H) φw () L2 μΓ () ⎛ ⎝⎞ ⎠2 k2ϕ1r () +ϕ2r () ~ x  w−x. (58) Setting v~ xin (48) and using (52) and (51), we obtain ~ α2φx () 2 L2 μΓ () ~ x2≤~ β2φx () L2 μΓ () uL2 μ(Γ;H)~ x. Thus, ~ x≤k2uL2 μ(Γ;H) φx () L2 μΓ () ≤k2uL2 μ(Γ;H) εk2 s. (59) It holds w w−~ x ~ x≤2w−~ x. Then, using (58) and (59),we deduce w w−~ x ~ x ≤2k2 s2ϕ1r () +ϕ2r () k2 s  w−x. That is, w−~ x ~ x ≤Δw−x,with Δgiven by (45). (60) Estimate (46) follows from this last inequality for x=w n ,assuming that w n ∈B H (w,r). Assuming w 0 ∈B H (w,r) this recursively proves that all the w n are in B H (w,r). Furthermore, suppose that there exists another solution to (35) with norm one in the ball B H (w,r), w*. In this case, estimating (60) for x=w*implies w−w*≤Δw−w* because ~ xT(w*)w*. Then w=w*, and there is uniqueness of solution with norm one in the ball B H (w,r). Finally, (47) follows from (46) by recurrence. Let us now consider method (41). In this case ~ x^ T(x),by (42), is the solution of the problem Γ φx, γ  2a~ x, v;γ  dμγ  Γ φx, γ  〈fγ ,v〉dμγ ,∀v∈H. (61) To estimate w−~ x ~ x, from problems (35) and (61), we obtain Γ φw, γ  2aw−~ x, v;γ  dμγ  −Γ φw, γ  2−φx, γ  2  a~ x, v;γ  dμγ  +Γ φw, γ  −φx, γ  〈fγ  ,v〉dμγ ,∀v∈H. (62) As 〈f(γ), v〉=a(u(γ), v;γ)≤βu(γ)v,∀v∈H, setting vw−~ x, we have TABLE 1 | Convergence rates of the PI algorithm for the first TN modes. n Mode i = 1 Mode i = 2 Mode i = 3 wn 1−wn 1−1rn 1wz 2n−wn−1 2rn 2wn 3−wn−1 3rn 3 1 7,4844E-01 —4,6117E-01 —3,8004E-02 — 2 9,6123E-02 7,78 8,3125E-02 5,54 5,9211E-02 0,641 3 4,1109E-03 23,38 8,3769E-03 9,92 6,8704E-03 8,61 4 1,7378E-04 23,65 7,9166E-04 10,58 7,5651E-04 9,08 5 7,3741E-06 23,56 7,4417E-05 10,63 8,2797E-04 9,13 6 3,1308E-07 23,55 6,9923E-06 10,64 9,0557E-06 9,14 7——6,5699E-07 10,64 9,9037E-07 9,14 Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696817 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric αφw () 2 L2 μΓ () w−~ x2≤Γ φw, γ  2aw−~ x, w −~ x;γ  dμγ ≤ βφ2w () −φ2x () L1 μΓ () ~ x+φw () −φx () L2 μ(Γ)uL2 μ(Γ;H)  w−~ x. (63) The functions φ(w) and φ(x) have the same expressions for methods (37) and (41). Moreover, setting v~ xin (61) αφw () 2 L2 μΓ () ~ x2≤Γ φw, γ  2a~ x, ~ x;γ  dμγ  Γ φw, γ  〈fγ  ,~ x〉dμγ  Γ auγ ,φw, γ  ~ x () ;γ  dμγ  ≤βuL2 μ(Γ;H)φw () L2 μΓ () ~ x. (64) Hence, ~ x≤β α s. Then, similarly to (58), we obtain w−~ x≤β α s2ϕ1r () +ϕ2r () ~ x  w−x ≤β α s2ϕ1r () +ϕ2r ()β α s  w−x. (65) As w−~ x ~ x≤2w−~ x, the conclusion follows as for method (37). Remark 2 The optimal convergence rate Δcorresponds to k =1 and λ=1, that is, ~ α~ βand α=β.Asαand βare predetermined, the optimal convergence rate can only be obtained with method (37). When the inner products (·,·) γ =(·,·) γ,+ and thus the operator A γ =A γ,+ , introduced in Section 3 are used to construct the optimal targeted subspaces, it satisfies, by (15), Aγ,+v2 γ,+v2 H,γ,∀v∈H. Then, from (44), choosing all the inner products · H,γ equal to the reference inner product ·,we obtain ~ α~ β1. Therefore,the convergence rate is optimal for method (37) with this choice. It can TABLE 2 | Convergence rates of the PI algorithm for the first STN modes. n Mode i = 1 Mode i = 2 Mode i = 3 wn 1−wn−1 1rn 1wn 2−wn−1 2rn 2wn 3−wn−1 3rn 3 1 8,0749E-1 —5,0427E-1 —3,9947E-1 — 2 1,4643E-1 5,51 1,7004E-1 2,96 1,1811E-1 3,38 3 9.8900E-3 14,08 4,0532E-2 4,19 8,9310E-2 1,32 4 6,0098E-4 16,45 8,8019E-3 4,60 4,1810E-2 2,13 5 3,6271E-5 16,56 1,8704E-3 4,70 1,8745E-2 2,23 6 2,1882E-6 16,57 3,9656E-4 4,72 8,2109E-3 2,28 7 1,3201E-7 16,57 8,3595E-5 4,73 3,5573E-3 2,30 8 7,9636E-9 16,57 1,7660E-5 4,73 1,5337E-3 2,31 9————6,5980E-4 2,32 FIGURE 1 | Solutions of problem (66) by the LRTD (called FLSTD within the figure) expansion computed with the TN method. Left γ= 2.7375 Right γ= 49.7624. FIGURE 2 | Convergence history of the PGD, TN, and STN series for Test 2. TABLE 3 | Numerical behavior of the PGD, STN, and TN methods for Test 2. Methods Errors L2 μ(Γ,L2(Ω)) Errors L2 μ(Γ,H1 0(Ω)) PGD u−u 94 = 9.69e−7u−u 94 = 3.30e−7 STN u−u 24 = 8.21e−7u−u 24 = 3.83e−7 TN u−u 12 = 6.31e−7u−u 12 = 3.52e−7 Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696818 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric be interpreted as preconditioning of the problem to solve, similar to classical preconditioning to accelerate convergence in solving linear systems (For example, see [24]). Remark 3 If we intend to compute a mode of order i ≥2,Theorem 6.1 and Remark 2also hold, replacing fby the residual fif−Aui−1 and uby the error u−u i−1 ,whereu i−1 is defined by (31). 7 NUMERICAL TESTS In this section, we discuss the numerical results obtained with the methods TN (38)–(39) and STN (41)–(42) to solve some non-symmetric second-order PDEs. Our purpose, on the one hand, is to confirm the theoretical results stated in Theorem 4.1 and Theorem 6.1, and on the other hand, to compare the practical performances of these methods with the standard PGD. We consider a parametric 1D advection–diffusion problem with fixed constant advection velocity β, γβu′−u″γfin Ω0,10 () u0 () 1,u′10 () 0 (66) We assume that β= 1, and then γis the Péclet number. The source term is f= 1/500. We have set Γ= [2.5, 50]; then, there is a large asymmetry of the advection–diffusion operator. Once the nonhomogenous boundary condition at x=0 is lifted, problem (66) is formulated under the general framework (1) when the space Hand the bilinear form a(·,·) are given by H={v∈H 1 (Ω)|v(0) = 0 }, and au,v;γ  γβu′,v  +u′,v′  . (67) We endow space Hwith the H1 0(Ω)norm (that we still denote ·), which is equivalent to the H 1 (Ω) norm on H. In practice, we replace the continuous problem (1) by an approximated one on a finite element space H h formed by piecewise affine elements. In addition, the integrals on Γare approximated by a quadrature formula constructed on a subdivision of Γinto Msubintervals, Γ ψγ dμγ ≃IMψ  M i1 ωiψγ i  . This is equivalent to approximating the Lebesgue measure μ by a discrete measure μ Δ located at the nodes of the discrete set Γ Δ ={γ i ,i=1,... ,M}withweightsω i ,i=1,... ,M. Consequently, all the theoretical results obtained in the previous sections apply, by replacing the L2 μ(Γ,H)space by L2 μΔ(ΓΔ,H h).Itistobenotedthat v2 L2 μΔ ΓΔ,Hh () IMv2  . (68) In our computations, we have used the midpoint quadrature formula with M= 100 equally spaced subintervals of Γto construct I M and constructed H h with 300 subintervals of Ωof the same length. Test 1: This first experiment is intended to check the theoretical results on the convergence rate of the PI algorithm, stated in Theorem 6.1, for the TN and STN methods: we consider optimal targeted subspaces, in the sense of the standard H1 0(Ω)norm. That is, using w, v () H,γw, v () H1 0Ω () w′,v′  ,∀w, v ∈H1Ω () to define the mappings A γ,+ by (14) and the form  b*by(22). For each mode w i of the LRTD expansion (31), we have estimated the numerical convergence rate of the PI algorithm by rn+1 iwn i−wn−1 i wn+1 i−wn i. (69) Tables 1 and 2show the norm of the difference between two consecutive approximations and the ratios rn i. We display the results for the first three modes. We observe that the PI method converges with a nearly constant rate for each mode, in agreement with Theorem 6.1. The convergence rate is larger for the TN method, also as expected from this theorem. It is also noted that the convergence rates are smaller for higherorder modes. In Figure 1, we present the comparison between the solution obtained by finite elements for γ= 2.7375 and γ= 49.7625 and the truncated series sum for the TN, the results for the STN are similar. Test 2: In this test, we compare the convergence rates of PGD, TN, and STN methods to obtain the LRTD expansion (31) for the problem (66). Figure 2 displays the errors of the truncated series with respect to the number of modes, in norm L2 μ(Γ,H 1 0(Ω)).A spectral convergence may be observed for the three expansions. We observe that the convergence of the TN expansion, in terms of the number of modes needed to achieve an error level, is much faster than the convergence of the STN expansion, while this one is faster than the PGD one. This is clarified if we consider the number of modes required to achieve an error smaller than a given level. We FIGURE 3 | Comparison of CPU times to compute the TN, STN, and PGD expansions. Frontiers in Physics | www.frontiersin.org May 2022 | Volume 10 | Article 8696819 Chacón Rebollo et al. Low-Rank Approximation Non-Symmetric