Low-Rank Approximations for Parametric Non-Symmetric Elliptic Problems

In this study, we obtained low-rank approximations for the solution of parametric non-symmetric 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.

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.

The present study is motivated by [5], where the authors present and analyze a generalization of the previous study [16] when operator A depends 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) may be found in the studies [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 non-symmetric problems.

Let $(\mathrm{\Gamma },\mathcal{B},\mu )$ be a measure space, where we assume that the measure μ is σ-finite. Let H be a separable Hilbert space endowed with the scalar product (⋅, ⋅) and associated norm ‖ ⋅‖, denote by H′ the dual space of H and by ⟨⋅, ⋅⟩ the duality pairing between H′ and H. We will consider the Lebesgue space $L_{\mu }^2(\mathrm{\Gamma })$ and the Bochner space $L_{\mu }^2(\mathrm{\Gamma };H)$ and its dual space $L_{\mu }^2(\mathrm{\Gamma };H^{\prime })$, denoting ⟨⋅, ⋅⟩ as the duality between $L_{\mu }^2(\mathrm{\Gamma };H)$ and $L_{\mu }^2(\mathrm{\Gamma };H^{\prime })$. We are interested in solving the parametric family of variational problems:

$\left\{\begin{array}{c}\text{\hspace{0.17em}Find\hspace{0.17em}}u:\mathrm{\Gamma }\to H\text{\hspace{0.17em}such\hspace{0.17em}that\hspace{0.17em}}\hfill \\ a\left(u\left(\gamma \right),v;\gamma \right)=⟨f\left(\gamma \right),v⟩,\forall v\in H,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma },\hfill \end{array}\right.(1)$

where $a(\cdot ,\cdot ;\gamma ):H\times H\to \mathbb{R}$ is a parameter-dependent, possibly non-symmetric, bilinear form and f(γ) ∈ H′ is a parameter-dependent 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,

$a\left(w,v;\gamma \right)\le \beta \Vert w\Vert \Vert v\Vert ,\forall w,v\in H,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma },(2)$

$\alpha \Vert w{\Vert }^2\le a\left(w,w;\gamma \right),\forall w\in H,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.(3)$

By the Lax–Milgram theorem, problem (1) admits a unique solution μ-a. e. γ ∈ Γ. To treat the measurability of u with respect to γ, let us consider the problem:

Let f be a function that belongs to $L_{\mu }^2(\mathrm{\Gamma };H^{\prime })$ such that f(γ) = f(γ) μ-a. e. γ ∈ Γ,

$\left\{\begin{array}{c}\text{\hspace{0.17em}Find\hspace{0.17em}}\mathbf{u}\in L_{\mu }^2(\mathrm{\Gamma };H)\text{\hspace{0.17em}such\hspace{0.17em}that\hspace{0.17em}}\hfill \\ \bar{a}\left(\mathbf{u},\mathbf{v}\right)=\langle \mathbf{f},\mathbf{v}\rangle ,\forall \mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };H),\hfill \end{array}\right.(4)$

where $\bar{a}:L_{\mu }^2(\mathrm{\Gamma };H)\times L_{\mu }^2(\mathrm{\Gamma };H)\to \mathbb{R}$ is defined by

$\bar{a}\left(\mathbf{w},\mathbf{v}\right)={\int }_{\mathrm{\Gamma }}a\left(\mathbf{w}\left(\gamma \right),\mathbf{v}(\gamma );\gamma \right)d\mu \left(\gamma \right).(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])

$\mathbf{u}\left(\gamma \right)=u\left(\gamma \right),\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.$

We shall consider a symmetrized reformulation of the formulation (4). Let us consider a family of inner products in H, ${\{{(\cdot ,\cdot )}_{\gamma }\}}_{\gamma \in \mathrm{\Gamma }}$ which generate the norms ‖ ⋅‖_γ uniformly equivalent to the ‖ ⋅‖ norm and there exist α_H > 0 and β_H > 0 such that,

${\alpha }_H\Vert v\Vert \le \Vert v{\Vert }_{\gamma }\le {\beta }_H\Vert v\Vert \text{\hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em}}\gamma \in \mathrm{\Gamma },(6)$

considering the associated scalar product in $L_{\mu }^2(\mathrm{\Gamma };H)$ and ∫_Γ(⋅,⋅)_γ dμ(γ). As $\bar{a}$ is continuous and coercive on $L_{\mu }^2(\mathrm{\Gamma };H)$, there exists a unique isomorphism in $L_{\mu }^2(\mathrm{\Gamma };H)$ that we denote $\mathcal{A}$, such that

${\int }_{\mathrm{\Gamma }}{\left(\left(\mathcal{A}\mathbf{w}\right)\left(\gamma \right),\mathbf{v}\left(\gamma \right)\right)}_{\gamma }d\mu \left(\gamma \right)=\bar{a}\left(\mathbf{w},\mathbf{v}\right),\forall \mathbf{w},\mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };H).(7)$

Let us define a new bilinear form $\bar{b}:L_{\mu }^2(\mathrm{\Gamma };H)\times L_{\mu }^2(\mathrm{\Gamma };H)\to \mathbb{R}$ by

$\bar{b}\left(\mathbf{w},\mathbf{v}\right)≔\bar{a}\left(\mathbf{w},\mathcal{A}\mathbf{v}\right).(8)$

It is to be noted that $\bar{b}$ is symmetric, as it can be written as

$\bar{b}\left(\mathbf{w},\mathbf{v}\right)={\int }_{\mathrm{\Gamma }}{\left(\left(\mathcal{A}\mathbf{w}\right)\left(\gamma \right),\left(\mathcal{A}\mathbf{v}\right)\left(\gamma \right)\right)}_{\gamma }d\mu \left(\gamma \right).(9)$

Thus, as a consequence of (2) and (3), the form $\bar{b}$ defines a scalar product in $L_{\mu }^2(\mathrm{\Gamma };H)$ and generates a norm equivalent to the usual norm is this space.

Now, problem (4) is equivalent to

$\left\{\begin{array}{c}\text{\hspace{0.17em}Find\hspace{0.17em}}\mathbf{u}\in L_{\mu }^2(\mathrm{\Gamma };H)\text{\hspace{0.17em}such\hspace{0.17em}that\hspace{0.17em}}\hfill \\ \bar{b}\left(\mathbf{u},\mathbf{v}\right)=\langle \mathbf{f},\mathcal{A}\mathbf{v}\rangle ,\forall \mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };H).\hfill \end{array}\right..(10)$

We shall use this formulation to build optimal approximations of problem (1) on subspaces of finite-dimension, when the form $\bar{a}(\cdot ,\cdot )$ is not symmetric. For a given integer k ≥ 1, we denote by S_k the Grassmannian variety of H formed by its subspaces of dimension smaller than or equal to k and consider the problem:

$\underset{Z\in {\text{S}}_k}{\min }\bar{b}\left(\mathbf{u}-{\mathbf{u}}_Z,\mathbf{u}-{\mathbf{u}}_Z\right),(11)$

where u is the solution of problem (4) and u_Z is its approximation given by the Galerkin approximation of problem (10) in $L_{\mu }^2(\mathrm{\Gamma };Z)$:

$\left\{\begin{array}{c}\text{Find\hspace{0.17em}}{\mathbf{u}}_Z\in L_{\mu }^2(\mathrm{\Gamma };Z)\text{\hspace{0.17em}such\hspace{0.17em}that}\hfill \\ \bar{b}\left({\mathbf{u}}_Z,\mathbf{v}\right)=\langle \mathbf{f},\mathcal{A}\mathbf{v}\rangle ,\forall \mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };Z).\hfill \end{array}\right.(12)$

Then, for any $k\in \mathbb{N}$, an optimal subspace of dimension smaller than or equal to k is the best subspace of the family S_k that minimizes the error between the solution u and its Galerkin approximation u_Z on this subspace, with respect to the mean quadratic norm generated by $\bar{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 ${\{{(\cdot ,\cdot )}_{H,\gamma }\}}_{\gamma \in \mathrm{\Gamma }}$ on H, which generate norms ‖ ⋅‖_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 $\bar{b}$ defined by (8) actually is

$\bar{b}\left(\mathbf{w},\mathbf{v}\right)={\int }_{\mathrm{\Gamma }}{\left(\mathbf{w}\left(\gamma \right),\mathbf{v}\left(\gamma \right)\right)}_{H,\gamma }d\mu \left(\gamma \right).$

In this way, the optimal subspaces are the solution of the problem

$\underset{Z\in {\text{S}}_k}{\min }{\int }_{\mathrm{\Gamma }}\Vert \mathbf{u}\left(\gamma \right)-{\mathbf{u}}_Z\left(\gamma \right){\Vert }_{H,\gamma }^{2}d\mu \left(\gamma \right).(13)$

Let us consider the operators A_γ,⋆: H↦H and the bilinear forms (⋅,⋅)_γ,⋆ on H × H defined by

$a\left(w,A_{\gamma ,\star }v;\gamma \right)={\left(w,v\right)}_{H,\gamma },\forall w,v\in H,(14)$

${\left(w,v\right)}_{\gamma ,\star }={\left(A_{\gamma ,\star }^{-1}w,A_{\gamma ,\star }^{-1}v\right)}_{H,\gamma },\forall w,v\in H.(15)$

It is to be noted that A_γ,⋆ is an isomorphism on H and consequently (⋅,⋅)_γ,⋆ is an inner product on H. Due to (2) and (3), the norms generated by these inner products are uniformly equivalent to the reference norm in H. Moreover, by (14) and (15),

${\left(A_{\gamma ,\star }w,v\right)}_{\gamma ,\star }={\left(w,A_{\gamma ,\star }^{-1}v\right)}_{H,\gamma }=a\left(w,v;\gamma \right).(16)$

Let us consider now the inner product in $L_{\mu }^2(\mathrm{\Gamma };H)$ given by ∫_Γ(⋅,⋅)_γ,⋆ and the isomorphism ${\mathcal{A}}_{\star }:L_{\mu }^2(\mathrm{\Gamma };H)\mapsto L_{\mu }^2(\mathrm{\Gamma };H)$ defined by

${\int }_{\mathrm{\Gamma }}{\left(\left({\mathcal{A}}_{\star }\mathbf{w}\right)\left(\gamma \right),\mathbf{v}\left(\gamma \right)\right)}_{\gamma ,\star }d\mu \left(\gamma \right)=\bar{a}\left(\mathbf{w},\mathbf{v}\right),\forall \mathbf{w},\mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };H).(17)$

Then, it holds.

Lemma 1 Let A_γ : H → H be the continuous linear operators defined by

${\left(A_{\gamma }w,v\right)}_{\gamma }=a\left(w,v;\gamma \right),\forall w,v\in H,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.(18)$

Then, it holds

$\left(\mathcal{A}\mathbf{w}\right)\left(\gamma \right)=A_{\gamma }\mathbf{w}\left(\gamma \right),\forall \mathbf{w}\in H,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.(19)$

This result follows from a standard argument using the separability of space H that we omit for brevity. Then, by Lemma 1, we have

$\left({\mathcal{A}}_{\star }\mathbf{w}\right)\left(\gamma \right)=A_{\gamma ,\star }\mathbf{w}\left(\gamma \right),\forall \mathbf{w}\in L_{\mu }^2(\mathrm{\Gamma };H),\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.(20)$

Let us denote by ${\bar{b}}_{\star }$ the bilinear form on $L_{\mu }^2(\mathrm{\Gamma };H)$ given by

${\bar{b}}_{\star }\left(\mathbf{w},\mathbf{v}\right)=\bar{a}\left(\mathbf{w},{\mathcal{A}}_{\star }\mathbf{v}\right),\forall \mathbf{w},\mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };H).(21)$

Then, by (20) and (14),

${\bar{b}}_{\star }\left(\mathbf{w},\mathbf{v}\right)={\int }_{\mathrm{\Gamma }}a\left(\mathbf{w}\left(\gamma \right),A_{\gamma ,\star }\mathbf{v}\left(\gamma \right);\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}{\left(\mathbf{w}\left(\gamma \right),\mathbf{v}\left(\gamma \right)\right)}_{H,\gamma }d\mu \left(\gamma \right).(22)$

As a consequence, the optimal subspaces obtained by the least-squares problem (11) when $\bar{b}={\bar{b}}_{\star }$ satisfy (13).

Remark 1 When the forms a (⋅, ⋅; γ) are symmetric, if we choose

${\left(w,v\right)}_{H,\gamma }=a\left(w,v;\gamma \right),\forall w,v\in H,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma },(23)$

then A_γ,⋆ defined by (14) is the identity operator in H. From (21), it follows that ${\bar{b}}_{\star }=\bar{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 $\le k$, obtained by deflation. For all N ≥ 1, we approximate

$\mathbf{u}\simeq {\mathbf{u}}_N≔\sum _{i=1}^N{\mathbf{s}}_i,\text{\hspace{0.17em}with\hspace{0.17em}}{\mathbf{s}}_i\in L_{\mu }^2(\mathrm{\Gamma };H),(24)$

computed by the following algorithm:

Initialization: let be u₀ = 0.

Iteration: assuming u_i−1, known for any i ≥ 1, set e_i−1≔u − u_i−1. Then,

${\mathbf{u}}_i={\mathbf{u}}_{i-1}+{\mathbf{s}}_i,\text{\hspace{0.17em}where\hspace{0.17em}}{\mathbf{s}}_i={\left({\mathbf{e}}_{i-1}\right)}_W(25)$

with ${({\mathbf{e}}_{i-1})}_W$ the approximation of e_i−1 given by the problem (12) on an optimal approximation subspace W solution of the problem (11),

$W≔\underset{Z\in S_k}{\mathrm{argmin}}\bar{b}\left({\mathbf{e}}_{i-1}-{\left({\mathbf{e}}_{i-1}\right)}_Z,{\mathbf{e}}_{i-1}-{\left({\mathbf{e}}_{i-1}\right)}_Z\right).(26)$

It is to be noted that this algorithm does not need to know the solution u of problem (4) since e_i−1 is defined in terms of the current residual ${\mathbf{f}}_i=\mathbf{f}-\mathcal{A}{\mathbf{u}}_{i-1}$ by

$\left\{\begin{array}{c}{\mathbf{e}}_{i-1}\in L_{\mu }^2(\mathrm{\Gamma };H)\text{\hspace{0.17em}such\hspace{0.17em}that\hspace{0.17em}}\hfill \\ \bar{a}\left({\mathbf{e}}_{i-1},\mathbf{v}\right)=\langle \mathbf{f},\mathbf{v}\rangle -\bar{a}\left({\mathbf{u}}_{i-1},\mathbf{v}\right)≔\langle {\mathbf{f}}_i,\mathbf{v}\rangle \forall \mathbf{v}\in L_{\mu }^2(\mathrm{\Gamma };H)\hfill \end{array}\right.$

The convergence of the u_N to u is stated in Theorem 5.3 of [5] as the form $\bar{b}$ is an inner product in $L_{\mu }^2(\mathrm{\Gamma };H)$, with the generated norm equivalent to the standard one.

Theorem 4.1 The truncated series u_N determined by the deflation algorithm (24)–(26) satisfying

$\underset{N\to \infty }{\lim }\bar{b}\left(\mathbf{u}-{\mathbf{u}}_N,\mathbf{u}-{\mathbf{u}}_N\right)=0.$

Consequently, u_N strongly converges in $L_{\mu }^2(\mathrm{\Gamma };H)$ to the solution u of 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

${\mathbf{u}}_Z(\gamma )=\phi \left(\gamma \right)w,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma },\text{\hspace{0.17em}for\hspace{0.17em}some\hspace{0.17em}}\phi \in L_{\mu }^2(\mathrm{\Gamma }),w\in H.$

Then, the problem (11) can be written as (see [5], Sect. 6)

$\underset{v\in H,\psi \in L_{\mu }^2\left(\mathrm{\Gamma }\right)}{\min }J\left(v,\psi \right),(27)$

where

$J\left(v,\psi \right)=\bar{b}\left(\mathbf{u}-\psi \otimes v,\mathbf{u}-\psi \otimes v\right).(28)$

Any solution of problem (27) has to verify the following conditions:

Proposition 1 If $(w,\phi )\in H\times L_{\mu }^2(\mathrm{\Gamma })$ is a solution of problem (27), then it is also a solution of the following coupled nonlinear problem:

$\begin{array}{ccc}\hfill \bar{a}\left(\phi \otimes w,\phi \otimes A_{\gamma }v\right)& =\hfill & \langle \mathbf{f},\phi \otimes A_{\gamma }v\rangle \forall v\in H,\hfill \end{array}(29)$

$\begin{array}{ccc}\hfill \bar{a}\left(\phi \otimes w,\psi \otimes A_{\gamma }w\right)& =\hfill & \langle \mathbf{f},\psi \otimes A_{\gamma }w\rangle \forall \psi \in L_{\mu }^2(\mathrm{\Gamma }).\hfill \end{array}(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\times L_{\mu }^2(\mathrm{\Gamma })\to \mathbb{R}$ and 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 J is 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 u given by

$\mathbf{u}\simeq {\mathbf{u}}_N≔\sum _{i=1}^N{\phi }_i\otimes w_i,(31)$

where the modes $(w_i,{\phi }_i)\in H\times L_{\mu }^2(\mathrm{\Gamma })$ are recursively computed as a solution of the problem $({\mathbf{f}}_i≔\mathbf{f}-\mathcal{A}{\mathbf{u}}_{i-1})$ :

$\left\{\begin{array}{cc}\bar{a}\left({\phi }_i\otimes w_i,{\phi }_i\otimes A_{\gamma }v\right)\hfill & =\langle {\mathbf{f}}_i,{\phi }_i\otimes A_{\gamma }v\rangle \forall v\in H,\hfill \\ \bar{a}\left({\phi }_i\otimes w_i,\psi \otimes A_{\gamma }{w}_i\right)\hfill & =\langle {\mathbf{f}}_i,\psi \otimes A_{\gamma }{w}_i\rangle \forall \psi \in L_{\mu }^2(\mathrm{\Gamma }).\hfill \end{array}\right.(32)$

If this problem is solved in such a way that its solution is also a solution of

$\underset{v\in H,\psi \in L_{\mu }^2\left(\mathrm{\Gamma }\right)}{\min }{J}_i\left(v,\psi \right),\text{with\hspace{0.17em}}{J}_i\left(v,\psi \right)=\bar{b}\left(\mathbf{u}-{\mathbf{u}}_{i-1}-\psi \otimes v,\mathbf{u}-{\mathbf{u}}_{i-1}-\psi \otimes v\right),(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.

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), which we assume to admit a nonzero solution. For simplicity, the notation f will 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

${\int }_{\mathrm{\Gamma }}\phi \left(\gamma \right)a\left(w,A_{\gamma }w;\gamma \right)\psi \left(\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}⟨\mathbf{f}\left(\gamma \right),A_{\gamma }w⟩\psi \left(\gamma \right)d\mu \left(\gamma \right),\forall \psi \in L_{\mu }^2(\mathrm{\Gamma }),$

and then

$\phi \left(\gamma \right)=\phi \left(w,\gamma \right)≔\frac{⟨\mathbf{f}\left(\gamma \right),A_{\gamma }w⟩}{a\left(w,A_{\gamma }w;\gamma \right)}\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.(34)$

Thus, problem (29)–(30) consists in

$\left\{\begin{array}{c}\text{\hspace{0.17em}Find\hspace{0.17em}}w\in H\text{\hspace{0.17em}such\hspace{0.17em}that}\hfill \\ {\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(w,A_{\gamma }v;\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}\phi \left(w,\gamma \right)⟨\mathbf{f}\left(\gamma \right),A_{\gamma }v⟩d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}\right.(35)$

with $\phi (w,\cdot )\in L_{\mu }^2(\mathrm{\Gamma })$ defined by (34). For simplicity, we shall denote by φ(w) the function φ(w, ⋅).

Let us define the operator T : H → H that transforms w ∈ H into T(w) ∈ H solution to the problem

$\begin{array}{c}{\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(T\left(w\right),A_{\gamma }v;\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}\phi \left(w,\gamma \right)⟨\mathbf{f}\left(\gamma \right),A_{\gamma }v⟩d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}(36)$

T is well-defined from (2) and (3), and a solution of (35) is a fixed point of this operator. Moreover, as A_γ is linear,

$\phi \left(\lambda w\right)={\lambda }^{-1}\phi \left(w\right)\text{\hspace{0.17em}and\hspace{0.17em}then\hspace{0.17em}}T\left(\lambda w\right)=\lambda T\left(w\right),\forall \lambda \in \mathbb{R}.$

Thus, if (w, φ) is a solution to (29)–(30), then (λ w, λ⁻¹φ) 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⁰ ∈ H such that φ⁰ = φ(w⁰) is not zero in $L_{\mu }^2(\mathrm{\Gamma })$.

Iteration: Knowing wⁿ ∈ H, the following is computed :

$\begin{array}{cc}\hfill \left.\text{i}\right)\hfill & {\tilde{w}}^{n+1}=T\left(w^n\right)\in H\hfill \\ \hfill \left.\text{ii}\right)\hfill & w^{n+1}=\frac{{\tilde{w}}^{n+1}}{\Vert {\tilde{w}}^{n+1}\Vert }\in H\hfill \\ \hfill \left.\text{iii}\right)\hfill & {\phi }^{n+1}=\phi \left(w^{n+1}\right)\in L_{\mu }^2(\mathrm{\Gamma }).\hfill \end{array}(37)$

The next result states that this iterative procedure is well-defined.

Lemma 2 It is assumed that for some nonzero w ∈ H, it holds that $\phi (w)\in L_{\mu }^2(\mathrm{\Gamma })$ is not zero. Then $\tilde{w}=T(w)$ is not zero in H and $\phi (\tilde{w})$ is not zero in $L_{\mu }^2(\mathrm{\Gamma })$.

Proof First, by reduction to the absurd, it is assumed that T(w) = 0. From (36), we have

${\int }_{\mathrm{\Gamma }}\phi \left(w,\gamma \right)\langle \mathbf{f}\left(\gamma \right),A_{\gamma }v\rangle d\mu \left(\gamma \right)=0,\forall v\in H.$

In particular, for v = w and taking into account (34), we deduce that

${\int }_{\mathrm{\Gamma }}\frac{|\langle \mathbf{f}\left(\gamma \right),A_{\gamma }w\rangle {|}^2}{a\left(w,A_{\gamma }w;\gamma \right)}d\mu \left(\gamma \right)=0.$

Then,

$\langle \mathbf{f}\left(\gamma \right),A_{\gamma }w\rangle =0,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma },$

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 $\phi (\tilde{w})=0$. From (34), we have

$\langle \mathbf{f}\left(\gamma \right),A_{\gamma }\tilde{w}\rangle =0,\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.$

Then, setting $v=\tilde{w}$ in (36) and using (19), we obtain

$\bar{b}\left(\phi \left(w\right)\otimes \tilde{w},\phi \left(w\right)\otimes \tilde{w}\right)={\int }_{\mathrm{\Gamma }}a\left(\phi \left(w,\gamma \right)\tilde{w},\phi \left(w,\gamma \right)A_{\gamma }\tilde{w}\right)d\mu \left(\gamma \right)=0.$

As $\bar{b}$ is a scalar product in $L_{\mu }^2(\mathrm{\Gamma };H)$, this implies that

$\Vert \phi \left(w\right)\otimes \tilde{w}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}=\Vert \phi \left(w\right){\Vert }_{L_{\mu }^2(\mathrm{\Gamma })}\Vert \tilde{w}\Vert =0.$

We have already proven that $\tilde{w}\ne 0$. So, φ(w) has to be equal to zero, in contradiction with the initial hypothesis. Thus, our assumption is false and $\phi (\tilde{w})$ is not zero.

This result proves that if w⁰ and φ⁰ 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, H is a space of large finite-dimension and the integral in Γ is approximated by some quadrature rule with nodes ${\{{\gamma }_i\}}_{i=1}^M$. The method requires the computation of $A_{{\gamma }_i}v$ for all the elements v on 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

$\phi \left(w,\gamma \right)=\frac{⟨\mathbf{f}\left(\gamma \right),A_{\gamma ,\star }w⟩}{\Vert w{\Vert }_{H,\gamma }^2}\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.(38)$

In addition, the problem (36) that defines the operator T simplifies to problem:

$\begin{array}{c}{\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^2{\left(T\left(w\right),v\right)}_{H,\gamma }d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}\phi \left(w,\gamma \right)⟨\mathbf{f}\left(\gamma \right),A_{\gamma ,\star }v⟩d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}(39)$

We shall refer to the method (38)–(39) as the TN (targeted-norm) 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

$\bar{a}\left({\phi }_i\otimes w_i,{\phi }_i\otimes A_{{\gamma }_0}v\right)=\langle {\mathbf{f}}_i,{\phi }_i\otimes A_{{\gamma }_0}v\rangle \forall v\in H,$

for some γ₀ ∈ Γ, then this equation is equivalent to

$\bar{a}\left({\phi }_i\otimes w_i,{\phi }_i\otimes v\right)=\langle {\mathbf{f}}_i,{\phi }_i\otimes v\rangle \forall v\in 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ⁿ ∈ H, the following is computed:

$\begin{array}{cc}\hfill \left.\text{i}\right)\hfill & {\tilde{w}}^{n+1}=\widehat{T}\left(w^n\right)\in H\hfill \\ \hfill \left.ii\right)\hfill & w^{n+1}=\frac{{\tilde{w}}^{n+1}}{\Vert {\tilde{w}}^{n+1}\Vert }\in H\hfill \\ \hfill \left.\text{iii}\right)\hfill & {\phi }^{n+1}=\phi \left(w^{n+1}\right)\in L_{\mu }^2(\mathrm{\Gamma }).\hfill \end{array}(41)$

where $\widehat{T}(w)$ is computed by

$\begin{array}{c}{\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(\widehat{T}\left(w\right),v;\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}\phi \left(w,\gamma \right)⟨\mathbf{f}\left(\gamma \right),v⟩d\mu \left(\gamma \right),\forall v\in H,\hfill \end{array},(42)$

where φ(w, γ) is defined 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

$\phi \left(\gamma \right)≔\frac{⟨\mathbf{f}\left(\gamma \right),w⟩}{a\left(w,w;\gamma \right)}\mu -\text{a.e.\hspace{0.17em}}\gamma \in \mathrm{\Gamma }.$

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 φⁿ remain in the exterior of some ball of positive radius, say ɛ > 0, of $L_{\mu }^2(\mathrm{\Gamma })$, that is,

$\Vert {\phi }^n{\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}\ge \epsilon ,n=0,1,2,\cdots .(43)$

This is a working hypothesis that makes sense whenever the mode that we intend to compute is not zero, considering that the wⁿ is normalized.

From the definition of A_γ and (6), it holds

$\tilde{\alpha }\Vert v\Vert \le \Vert A_{\gamma }v{\Vert }_{\gamma }\le \tilde{\beta }\Vert v\Vert \text{\hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em}}\gamma \in \mathrm{\Gamma },\text{with\hspace{0.17em}}\tilde{\alpha }=\frac{\alpha }{{\beta }_H},\tilde{\beta }=\frac{\beta }{{\alpha }_H},(44)$

where α and β are given by (2) and (3), respectively, and α_H and β_H are defined in (6).

Let us define the function:

$\delta \left(r,s\right)=2\lambda k^2\left(1+\frac{2-r}{1-r}{k}^2\right)\left(1+\frac{2-r}{1-r}\lambda k^{2}s\right)s^2,$

$\text{where\hspace{0.17em}}k=\frac{\tilde{\beta }}{\tilde{\alpha }},\text{and}\lambda =\left\{\begin{array}{cc}{k}^2\hfill & \hfill \text{\hspace{0.17em}for\hspace{0.17em}method\hspace{0.17em}}\left(\text{37}\right)\hfill \\ \frac{\beta }{\alpha }\hfill & \hfill \text{\hspace{0.17em}for\hspace{0.17em}method\hspace{0.17em}}\left(\text{41}\right)\hfill \end{array}\right.$

It holds.

Theorem 6.1 It is assumed that (43) holds, and

$\mathrm{\Delta }=\delta \left(r,\bar{s}\right)<1,\mathrm{for\hspace{0.17em}some\hspace{0.17em}}r\in \left(\mathrm{0,1}\right)\mathrm{\hspace{0.17em}with\hspace{0.17em}}\bar{s}=\frac{\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}}{\epsilon }.(45)$

Then, there exists a unique solution with norm 1 w of problem (35); the sequence ${\{w^n\}}_{n\ge 1}$ computed by either method (37) or method (41) is contained in the ball B_H(w, r) if w⁰ ∈ B_H(w, r), and

$\Vert w-w^{n+1}\Vert \le \mathrm{\Delta }\Vert w-w^n\Vert ,\forall n\ge 0.(46)$

As a consequence, the sequence ${\{w^n\}}_{n\ge 1}$ that is defined by either method (37) or method (41) is strongly convergent to w with linear rate and the following error estimate holds:

$\Vert w-w^n\Vert \le {\mathrm{\Delta }}^n\Vert w-w^0\Vert ,\forall n\ge 1,\mathrm{whenever\hspace{0.17em}}{w}^0\in B_H\left(w,r\right).(47)$

Proof. Let us consider at first the method (37). Let x ∈ B_H(w, r) such that $\Vert \phi (x){\Vert }_{L_{\mu }^2(\mathrm{\Gamma })}\ge \epsilon$. Denote $\tilde{x}=T(x)$ which by the definition of operator T in (36) is the solution to the problem

$\begin{array}{c}{\int }_{\mathrm{\Gamma }}\phi {\left(x,\gamma \right)}^{2}a\left(\tilde{x},A_{\gamma }v;\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}\phi \left(x,\gamma \right)⟨\mathbf{f}\left(\gamma \right),A_{\gamma }v⟩d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}(48)$

We aim to estimate $\Vert w-\frac{\tilde{x}}{\Vert \tilde{x}\Vert }\Vert$. To do that, from problems (35) and (48), we obtain

$\begin{array}{ll}\hfill {\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(w-\tilde{x},A_{\gamma }v;\gamma \right)d\mu \left(\gamma \right)& =-{\int }_{\mathrm{\Gamma }}\left(\phi {\left(w,\gamma \right)}^2-\phi {\left(x,\gamma \right)}^2\right)a\left(\tilde{x},A_{\gamma }v;\gamma \right)d\mu \left(\gamma \right)\hfill \\ \hfill & +{\int }_{\mathrm{\Gamma }}\left(\phi \left(w,\gamma \right)-\phi \left(x,\gamma \right)\right)\hfill \\ \hfill & \langle \mathbf{f}\left(\gamma \right),A_{\gamma }v\rangle d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}(49)$

It holds

$a\left(y,A_{\gamma }z;\gamma \right)={\left(A_{\gamma }y,A_{\gamma }z\right)}_{\gamma }\le {\tilde{\beta }}^2\Vert y\Vert \Vert z\Vert ,\forall y,z\in H,(50)$

using (44). Thus,

$\langle \mathbf{f}\left(\gamma \right),A_{\gamma }v\rangle =a\left(u\left(\gamma \right),A_{\gamma }v;\gamma \right)\le {\tilde{\beta }}^2\Vert u\left(\gamma \right)\Vert \Vert v\Vert ,\forall v\in H.(51)$

Moreover,

$a\left(y,A_{\gamma }y;\gamma \right)={\left(A_{\gamma }y,A_{\gamma }y\right)}_{\gamma }\ge {\tilde{\alpha }}^2\Vert y{\Vert }^2\forall y\in H.(52)$

Setting $v=w-\tilde{x}$ in (62) and using (50)-(52), we have

$\begin{array}{c}{\tilde{\alpha }}^2\Vert \phi \left(w\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}^2\Vert w-\tilde{x}{\Vert }^2\le \hfill \\ {\tilde{\beta }}^2\left(\Vert {\phi }^2\left(w\right)-{\phi }^2\left(x\right){\Vert }_{L_{\mu }^1\left(\mathrm{\Gamma }\right)}\Vert \tilde{x}\Vert +\Vert \phi \left(w\right)-\phi \left(x\right){\Vert }_{L_{\mu }^2(\mathrm{\Gamma })}\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}\right)\Vert w-\tilde{x}\Vert .\hfill \end{array}(53)$

To bound the second term in the r. h s. of (53), from (34), it holds

$\begin{array}{ll}\hfill \phi \left(w,\gamma \right)-\phi \left(x,\gamma \right)& =\frac{a\left(u\left(\gamma \right),A_{\gamma }\left(w-x\right);\gamma \right)}{a\left(w,A_{\gamma }w;\gamma \right)}\hfill \\ \hfill & +\frac{a\left(x,A_{\gamma }\left(x-w\right);\gamma \right)+a\left(x-w,A_{\gamma }w;\gamma \right)}{a\left(w,A_{\gamma }w;\gamma \right)a\left(x,A_{\gamma }x;\gamma \right)}\hfill \\ \hfill & a\left(u\left(\gamma \right),A_{\gamma }x;\gamma \right).\hfill \end{array}$

Then, using (50) and (52),

$\begin{array}{ccc}|\phi \left(w,\gamma \right)-\phi \left(x,\gamma \right)|\hfill & \hfill \le \hfill & k^2\left(1+k^2\left(\frac{1}{\Vert x\Vert }+1\right)\right)\Vert u\left(\gamma \right)\Vert \Vert w-x\Vert \hfill \\ \hfill & \hfill \le \hfill & k^2\left(1+\frac{2-r}{1-r}{k}^2\right)\Vert u\left(\gamma \right)\Vert \Vert w-x\Vert ,\hfill \end{array}(54)$

where $k=\frac{\tilde{\beta }}{\tilde{\alpha }}$. In the last estimate, we have used that as x ∈ B_H(w, r),

‖w‖ − r ≤ ‖x‖ and then $\frac{1}{\Vert x\Vert }\le \frac{1}{1-r}$. Therefore,

$\Vert \phi \left(w\right)-\phi \left(x\right){\Vert }_{L_{\mu }^2(\mathrm{\Gamma })}\le {\varphi }_1\left(r\right)\Vert u{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}\Vert w-x\Vert ,(55)$

where

${\varphi }_1\left(r\right)=k^2\left(1+\frac{2-r}{1-r}{k}^2\right).$

It is to be noted that from (34),

$|\phi \left(z,\gamma \right)|=\left|\frac{a\left(u\left(\gamma \right),A_{\gamma }z;\gamma \right)}{a\left(z,A_{\gamma }z;\gamma \right)}\right|\le k^2\frac{\Vert u\left(\gamma \right)\Vert }{\Vert z\Vert },\forall z\in H.(56)$

Then, from (54) and (56)

$\begin{array}{ccc}|\phi {\left(w,\gamma \right)}^2-\phi {\left(x,\gamma \right)}^2|\hfill & \hfill \le \hfill & |\phi \left(w,\gamma \right)-\phi \left(x,\gamma \right)||\phi \left(w,\gamma \right)+\phi \left(x,\gamma \right)|\hfill \\ \hfill & \hfill \le \hfill & k^2{\varphi }_1\left(r\right)\left(1+\frac{1}{\Vert x\Vert }\right)\Vert u\left(\gamma \right){\Vert }^2\Vert w-x\Vert .\hfill \end{array}$

Hence,

$\Vert \phi {\left(w\right)}^2-\phi {\left(x\right)}^2{\Vert }_{L_{\mu }^1(\mathrm{\Gamma })}\le {\varphi }_2\left(r\right)\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}^2\Vert w-x\Vert ,(57)$

where

${\varphi }_2\left(r\right)=k^2{\varphi }_1\left(r\right)\frac{2-r}{1-r}.$

Combining (53) with (55) and (57), we deduce

$\Vert w-\tilde{x}\Vert \le {\left(\frac{\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}}{\Vert \phi \left(w\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}}\right)}^2{k}^2\left({\varphi }_1\left(r\right)+{\varphi }_2\left(r\right)\Vert \tilde{x}\Vert \right)\Vert w-x\Vert .(58)$

Setting $v=\tilde{x}$ in (48) and using (52) and (51), we obtain

${\tilde{\alpha }}^2\Vert \phi \left(x\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}^2\Vert \tilde{x}{\Vert }^2\le {\tilde{\beta }}^2\Vert \phi \left(x\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}\Vert \tilde{x}\Vert .$

Thus,

$\Vert \tilde{x}\Vert \le k^2\frac{\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}}{\Vert \phi \left(x\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}}\le k^2\frac{\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}}{\epsilon }=k^2\bar{s}.(59)$

It holds $‖\frac{w}{\Vert w\Vert }-\frac{\tilde{x}}{\Vert \tilde{x}\Vert }‖\le 2\Vert w-\tilde{x}\Vert .$ Then, using (58) and (59), we deduce

$‖\frac{w}{\Vert w\Vert }-\frac{\tilde{x}}{\Vert \tilde{x}\Vert }‖\le 2k^2{\bar{s}}^2\left({\varphi }_1\left(r\right)+{\varphi }_2\left(r\right)k^2\bar{s}\right)\Vert w-x\Vert .$

That is,

$‖w-\frac{\tilde{x}}{\Vert \tilde{x}\Vert }‖\le \mathrm{\Delta }\Vert w-x\Vert ,\text{with\hspace{0.17em}}\mathrm{\Delta }\text{given\hspace{0.17em}by\hspace{0.17em}}\left(\mathbf{\text{45}}\right).(60)$

Estimate (46) follows from this last inequality for x = wⁿ, assuming that wⁿ ∈ B_H(w, r). Assuming w⁰ ∈ B_H(w, r) this recursively proves that all the wⁿ 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

$\Vert w-w^{*}\Vert \le \mathrm{\Delta }\Vert w-w^{*}\Vert$

because $\tilde{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 $\tilde{x}=\widehat{T}(x)$, by (42), is the solution of the problem

$\begin{array}{c}{\int }_{\mathrm{\Gamma }}\phi {\left(x,\gamma \right)}^{2}a\left(\tilde{x},v;\gamma \right)d\mu \left(\gamma \right)={\int }_{\mathrm{\Gamma }}\phi \left(x,\gamma \right)⟨\mathbf{f}\left(\gamma \right),v⟩d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}(61)$

To estimate $‖w-\frac{\tilde{x}}{\Vert \tilde{x}\Vert }‖$, from problems (35) and (61), we obtain

$\begin{array}{c}{\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(w-\tilde{x},v;\gamma \right)d\mu \left(\gamma \right)=\hfill \\ -{\int }_{\mathrm{\Gamma }}\left(\phi {\left(w,\gamma \right)}^2-\phi {\left(x,\gamma \right)}^2\right)a\left(\tilde{x},v;\gamma \right)d\mu \left(\gamma \right)\hfill \\ +{\int }_{\mathrm{\Gamma }}\left(\phi \left(w,\gamma \right)-\phi \left(x,\gamma \right)\right)\langle \mathbf{f}\left(\gamma \right),v\rangle d\mu \left(\gamma \right),\forall v\in H.\hfill \end{array}(62)$

As ⟨f(γ), v⟩ = a (u(γ), v; γ) ≤ β ‖u(γ)‖ ‖v‖, ∀ v ∈ H, setting $v=w-\tilde{x}$, we have

$\begin{array}{ccc}\hfill & \alpha \Vert \phi \left(w\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}^2\Vert w-\tilde{x}{\Vert }^2\le {\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(w-\tilde{x},w-\tilde{x};\gamma \right)d\mu \left(\gamma \right)\le \hfill & \hfill \\ \hfill & \beta \left(\Vert {\phi }^2\left(w\right)-{\phi }^2\left(x\right){\Vert }_{L_{\mu }^1\left(\mathrm{\Gamma }\right)}\Vert \tilde{x}\Vert +\Vert \phi \left(w\right)-\phi \left(x\right){\Vert }_{L_{\mu }^2(\mathrm{\Gamma })}\Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}\right)\Vert w-\tilde{x}\Vert .\hfill & \hfill \end{array}(63)$

The functions φ(w) and φ(x) have the same expressions for methods (37) and (41).

Moreover, setting $v=\tilde{x}$ in (61)

$\begin{array}{ll}\hfill \alpha \Vert \phi \left(w\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}^2\Vert \tilde{x}{\Vert }^2& \le {\int }_{\mathrm{\Gamma }}\phi {\left(w,\gamma \right)}^{2}a\left(\tilde{x},\tilde{x};\gamma \right)d\mu \left(\gamma \right)\hfill \\ \hfill & ={\int }_{\mathrm{\Gamma }}\phi \left(w,\gamma \right)⟨f\left(\gamma \right),\tilde{x}⟩d\mu \left(\gamma \right)\hfill \\ \hfill & ={\int }_{\mathrm{\Gamma }}a\left(\mathbf{u}\left(\gamma \right),\phi \left(w,\gamma \right)\left(\tilde{x}\right);\gamma \right)d\mu \left(\gamma \right)\hfill \\ \hfill & \le \beta \Vert \mathbf{u}{\Vert }_{L_{\mu }^2(\mathrm{\Gamma };H)}\Vert \phi \left(w\right){\Vert }_{L_{\mu }^2\left(\mathrm{\Gamma }\right)}\Vert \tilde{x}\Vert .\hfill \end{array}(64)$

Hence, $\Vert \tilde{x}\Vert \le \frac{\beta }{\alpha }\bar{s}.$ Then, similarly to (58), we obtain

$\begin{array}{ccc}\hfill \Vert w-\tilde{x}\Vert & \le \hfill & \frac{\beta }{\alpha }{\bar{s}}^2\left({\varphi }_1\left(r\right)+{\varphi }_2\left(r\right)\Vert \tilde{x}\Vert \right)\Vert w-x\Vert \hfill \\ \hfill & \le \hfill & \frac{\beta }{\alpha }{\bar{s}}^2\left({\varphi }_1\left(r\right)+{\varphi }_2\left(r\right)\frac{\beta }{\alpha }\bar{s}\right)\Vert w-x\Vert .\hfill \end{array}(65)$

As $‖w-\frac{\tilde{x}}{\Vert \tilde{x}\Vert }‖\le 2‖w-\tilde{x}‖$, the conclusion follows as for method (37).

Remark 2 The optimal convergence rate Δ corresponds to k = 1 and λ = 1, that is, $\tilde{\alpha }=\tilde{\beta }$ 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),

$\Vert A_{\gamma ,\star }v{\Vert }_{\gamma ,\star }^2=\Vert v{\Vert }_{H,\gamma }^2,\forall v\in H.$

Then, from (44), choosing all the inner products ‖ ⋅‖_H,γ equal to the reference inner product ‖ ⋅‖, we obtain $\tilde{\alpha }=\tilde{\beta }=1$. Therefore, the convergence rate is optimal for method (37) with this choice. It can 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 2 also hold, replacing f by the residual ${\mathbf{f}}_i=\mathbf{f}-\mathcal{A}{\mathbf{u}}_{i-1}$ and u by the error u − u_i−1, where u_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 β,

$\left\{\begin{array}{ccc}\gamma \beta u^{\prime }-u^{″}=\gamma f\hfill & \hfill \text{in}\hfill & \mathrm{\Omega }=\left(\mathrm{0,10}\right)\hfill \\ u\left(0\right)=1,u^{\prime }\left(10\right)=0\hfill & \hfill \hfill & \hfill \end{array}\right.(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 H and the bilinear form a (⋅, ⋅) are given by H = {v ∈ H¹(Ω) |v (0) = 0 }, and

$a\left(u,v;\gamma \right)=\gamma \left(\beta u^{\prime },v\right)+\left(u^{\prime },v^{\prime }\right).(67)$

We endow space H with the $H_0^1(\mathrm{\Omega })$ norm (that we still denote ‖ ⋅‖), which is equivalent to the H¹(Ω) 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 M subintervals,

${\int }_{\mathrm{\Gamma }}\psi \left(\gamma \right)d\mu \left(\gamma \right)\simeq I_M\left(\psi \right)=\sum _{i=1}^M{\omega }_i\psi \left({\gamma }_i\right).$

This is equivalent to approximating the Lebesgue measure μ by a discrete measure μ_Δ located at the nodes of the discrete set Γ_Δ = {γ_i, i = 1, … , M} with weights ω_i, i = 1, … , M. Consequently, all the theoretical results obtained in the previous sections apply, by replacing the $L_{\mu }^2(\mathrm{\Gamma },H)$ space by $L_{{\mu }_{\mathrm{\Delta }}}^2({\mathrm{\Gamma }}_{\mathrm{\Delta }},H_h)$. It is to be noted that

$\Vert \mathbf{v}{\Vert }_{L_{{\mu }_{\mathrm{\Delta }}}^2\left({\mathrm{\Gamma }}_{\mathrm{\Delta }},H_h\right)}^2=I_M\left(\Vert \mathbf{v}{\Vert }^2\right).(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 $H_0^1(\mathrm{\Omega })$ norm. That is, using

${\left(w,v\right)}_{H,\gamma }={\left(w,v\right)}_{H_0^1\left(\mathrm{\Omega }\right)}=\left(w^{\prime },v^{\prime }\right),\forall w,v\in H^1\left(\mathrm{\Omega }\right)$

to define the mappings A_γ,⋆ by (14) and the form ${\bar{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

$r_i^{n+1}=\frac{\Vert w_i^n-w_i^{n-1}\Vert }{\Vert w_i^{n+1}-w_i^n\Vert }.(69)$

Tables 1 and 2 show the norm of the difference between two consecutive approximations and the ratios $r_i^n$. We display the results for the first three modes.

TABLE 1

TABLE 1. Convergence rates of the PI algorithm for the first TN modes.

TABLE 2

TABLE 2. Convergence rates of the PI algorithm for the first STN 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 higher-order 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.

FIGURE 1

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.

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 $L_{\mu }^2(\mathrm{\Gamma },H_0^1(\mathrm{\Omega }))$. 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 display these numbers for an error level of 10^–6 in Table 3, where much more modes are needed by the PGD expansion. The TN and STN methods, thus, appear to be well-adapted to fit the asymmetry of the operator.

FIGURE 2

FIGURE 2. Convergence history of the PGD, TN, and STN series for Test 2.

TABLE 3

TABLE 3. Numerical behavior of the PGD, STN, and TN methods for Test 2.

Finally, we compare the CPU times required by the three methods. By construction, it is clear that to compute every single iteration of the PI method, the TN method is much more expensive since it involves the calculation of the A_γ,* operator for each finite element base function. However, due to the fast convergence of the associated LRTD expansion, it is less expensive than PGD to compute the expansion. Figure 3 displays the CPU times for the TN, STN, and PGD methods as a function of the number of subintervals M considered in the partition of Γ. The STN method is more expensive than the PGD method; this arises due to the small convergence rate of the PI algorithm with the STN method. However, the TN method is less expensive than the PGD one, requiring approximately half the CPU time.

FIGURE 3

FIGURE 3. Comparison of CPU times to compute the TN, STN, and PGD expansions.

8 Conclusion

In this study, we have proposed a new low-rank tensorized decomposition (LRTD) to approximate the solution of parametric non-symmetric elliptic problems, based on symmetrization of the problem.

Each mode of the series is characterized as a solution to a calculus of variation problem that yields an optimal finite-dimensional subspace, in the sense that it minimizes the error between the parametric solution and its approximation on this subspace, with respect to a preset mean quadratic norm. We have proven that the truncated expansion given by the deflation algorithm strongly converges to the parametric solution in the mean quadratic norm.

The minimization problems to compute the rank-one optimal modes are solved by the power iteration algorithm. We have proven that this method is locally linearly convergent when the initial data are close enough to an optimal mode. We also have identified an optimal symmetrization that provides the best convergence rates of the PI algorithm, with respect to the preset mean quadratic norm.

Furthermore, we have presented 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 compare advantageously to those required by the PGD expansion.

In this study, we have focused on rank-one tensorized decompositions. In our forthcoming research, we intend to extend the analysis to ranks k ≥ 2. This requires solving minimization problems on a Grassmann variety to compute the LRTD modes. We will also work on the solution of higher-dimensional non-symmetric elliptic problems by the method introduced in order to reduce the computation times as these increase with the dimension of the approximation spaces.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Author Contributions

TC and IS performed the theoretical derivations while MG did the computations.

Funding

The research of TC and IS has been partially funded by PROGRAMA OPERATIVO FEDER ANDALUCIA 2014-2020 grant US-1254587, which in addition has covered the cost of this publication. The research of MG is partially funded by the Spanish Research Agency - EU FEDER Fund grant RTI2018-093521-B-C31.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Ammar A, Chinesta F, Falcó A. On the Convergence of a Greedy Rank-One Update Algorithm for a Class of Linear Systems. Arch Computat Methods Eng (2010) 17:473–86. Number 4. doi:10.1007/s11831-010-9048-z

CrossRef Full Text | Google Scholar

2. Ammar A, Mokdad B, Chinesta F, Keunings R. A New Family of Solvers for Some Classes of Multidimensional Partial Differential Equations Encountered in Kinetic Theory Modeling of Complex Fluids. J Non-Newtonian Fluid Mech (2006) 139:153–76. doi:10.1016/j.jnnfm.2006.07.007

CrossRef Full Text | Google Scholar

3. Azaïez M, Chacón Rebollo T, Gómez Mármol M. On the Computation of Proper Generalized Decomposition Modes of Parametric Elliptic Problems. SeMA (2020) 77:59–72. doi:10.1007/s40324-019-00198-7

CrossRef Full Text | Google Scholar

4. Azaïez M, Belgacem FB, Rebollo TC. Error Bounds for POD Expansions of Parameterized Transient Temperatures. Comp Methods Appl Mech Eng (2016) 305:501–11. doi:10.1016/j.cma.2016.02.016

CrossRef Full Text | Google Scholar

5. Azaïez M, Belgacem FB, Casado-Díaz J, Rebollo TC, Murat F. A New Algorithm of Proper Generalized Decomposition for Parametric Symmetric Elliptic Problems. SIAM J Math Anal (2018) 50(5):5426–45. doi:10.1137/17m1137164

CrossRef Full Text | Google Scholar

6. ten Berge JMF, de Leeuw J, Kroonenberg PM. Some Additional Results on Principal Components Analysis of Three-Mode Data by Means of Alternating Least Squares Algorithms. Psychometrika (1987) 52:183–91. doi:10.1007/bf02294233

CrossRef Full Text | Google Scholar

7. Bulut H, Akkilic AN, Khalid BJ. Soliton Solutions of Hirota Equation and Hirota-Maccari System by the (m+1/G’)-expansion Method. Adv Math Models Appl (2021) 6(1):22–30.

Google Scholar

8. Le Bris C, Lelièvre T, Maday Y. Results and Questions on a Nonlinear Approximation Approach for Solving High-Dimensional Partial Differential Equations. Constr Approx (2009) 30(3):621–51. doi:10.1007/s00365-009-9071-1

CrossRef Full Text | Google Scholar

9. Cancès E, Lelièvre T, Ehrlacher V. Convergence of a Greedy Algorithm for High-Dimensional Convex Nonlinear Problems. Math Models Methods Appl Sci (2011) 21(12):2433–67. doi:10.1142/s0218202511005799

CrossRef Full Text | Google Scholar

10. Cancès E, Lelièvre T, Ehrlacher V. Greedy Algorithms for High-Dimensional Non-symmetric Linear Problems. Esaim: Proc (2013) 41:95–131. doi:10.1051/proc/201341005

CrossRef Full Text | Google Scholar

11. Chinesta F, Ammar A, Cueto E. Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models. Arch Computat Methods Eng (2010) 17(4):327–50. doi:10.1007/s11831-010-9049-y

CrossRef Full Text | Google Scholar

12. Espig M, Hackbusch W, Rohwedder T, Schneider R. Variational Calculus with Sums of Elementary Tensors of Fixed Rank. Numer Math (2012) 122:469–88. doi:10.1007/s00211-012-0464-x

CrossRef Full Text | Google Scholar

13. Espig M, Hackbusch W. A Regularized Newton Method for the Efficient Approximation of Tensors Represented in the Canonical Tensor Format. Numer Math (2012) 122:489–525. doi:10.1007/s00211-012-0465-9

CrossRef Full Text | Google Scholar

14. Espig M, Hackbusch W, Khachatryan A. On the Convergence of Alternating Least Squares Optimisation in Tensor Format Representations. arXiv:1506.00062v1 [math.NA] (2015).

Google Scholar

15. Espig M, Hackbusch W, Litvinenko A, Matthies HG, Zander E. Iterative Algorithms for the post-processing of High-Dimensional Data. J Comput Phys (2020) 410:109–396. doi:10.1016/j.jcp.2020.109396

CrossRef Full Text | Google Scholar

16. Falcó A, Nouy A. A Proper Generalized Decomposition for the Solution of Elliptic Problems in Abstract Form by Using a Functional Eckart-Young Approach. J Math Anal Appl (2011) 376:469–80. doi:10.1016/j.jmaa.2010.12.003

CrossRef Full Text | Google Scholar

17. Falcó A, Nouy A. Proper Generalized Decomposition for Nonlinear Convex Problems in Tensor Banach Spaces. Numer Math (2012) 121:503–30. doi:10.1007/s00211-011-0437-5

CrossRef Full Text | Google Scholar

18. Nouy A. Low-rank Tensor Methods for Model Order Reduction. In: Handbook of Uncertainty Quantification (Roger GhanemDavid HigdonHouman Owhadi. Eds. Philadelphia, PA: Springer (2017). doi:10.1007/978-3-319-12385-1_21

CrossRef Full Text | Google Scholar

19. Nouy A. Low-rank Methods for High-Dimensional Approximation and Model Order Reduction. In: P Benner, A Cohen, M Ohlberger, and K Willcox, editors. Model Reduction and Approximations. Philadelphia, PA: SIAM (2017).

Google Scholar

20. Kiers HAL. An Alternating Least Squares Algorithm for PARAFAC2 and Three-Way DEDICOM. Comput Stat Data Anal (1993) 16:103–18. doi:10.1016/0167-9473(93)90247-q

CrossRef Full Text | Google Scholar

21. Kroonenberg PM, de Leeuw J. Principal Component Analysis of Three-Mode Data of using Alternating Least Squares Algorithms. Psychometrika (1980) 45:69–97. doi:10.1007/bf02293599

CrossRef Full Text | Google Scholar

22. Ladévèze P. Nonlinear Computational Structural Mechanics. New Approaches and Non-incremental Methods of Calculation. Berlin: Springer (1999).

Google Scholar

23. Mohlenkamp MJ. Musings on Multilinear Fitting. Linear algebra and its applications (2013) 438(2): 834–52. doi:10.1016/j.laa.2011.04.019

CrossRef Full Text | Google Scholar

24. Rasheed SM, Nachaoui A, Hama MF, Jabbar AK. Regularized and Preconditioned Conjugate Gradient Like-Methods Methods for Polynomial Approximation of an Inverse Cauchy Problem. Adv Math Models Appl (2021) 6(2):89–105.

Google Scholar

25. Uschmajew A. Local Convergence of the Alternating Least Squares Algorithm for Canonical Tensor Approximation. SIAM J Matrix Anal Appl (2012) 33(2):639–52. doi:10.1137/110843587

CrossRef Full Text | Google Scholar

26. Wang L, Chu MT. On the Global Convergence of the Alternating Least Squares Method for Rank-One Approximation to Generic Tensors. SIAM J Matrix Anal Appl (2012) 35(3):1058–72.

Google Scholar

27. Zhang T, Golub GH. Rank-one Approximation to High Order Tensors. SIAM J Matrix Anal Appl (2001) 23(2):534–50. doi:10.1137/s0895479899352045

CrossRef Full Text | Google Scholar