Thomas Boiveau

Virginie Ehrlacher

Alexandre Ern

Anthony Nouy

Low-rank approximation of linear parabolic equations by space-time tensor Galerkin methods *

Keywords: Parabolic equations, Tensor methods, Proper Generalized Decomposition, Greedy algorithm. AMS. 65M12, 65M22, 35K20

We devise a space-time tensor method for the low-rank approximation of linear parabolic evolution equations. The proposed method is a stable Galerkin method, uniformly in the discretization parameters, based on a Minimal Residual formulation of the evolution problem in Hilbert-Bochner spaces. The discrete solution is sought in a linear trial space composed of tensors of discrete functions in space and in time and is characterized as the unique minimizer of a discrete functional where the dual norm of the residual is evaluated in a space semi-discrete test space. The resulting global space-time linear system is solved iteratively by a greedy algorithm. Numerical results are presented to illustrate the performance of the proposed method on test cases including non-selfadjoint and time-dependent differential operators in space. The results are also compared to those obtained using a fully discrete Petrov-Galerkin setting to evaluate the dual residual norm.

Introduction

The goal of this work is to devise a space-time tensor method for the low-rank approximation of the solution to linear parabolic evolution equations. The method we propose has two salient features. First, it is a stable Galerkin method, uniformly with respect to the discretization parameters in space and in time, leading to quasi-optimal error estimates in the natural norms of the problem as specified below. Second, the method is global in time (and in space), and the approximate solution is iteratively constructed by solving alternatively global problems 1 in space and in time. More precisely, at iteration m ∈ N * , the space-time approximate solution u m (x, t) is of the form

u m (x, t) = 1≤n≤m v n (x)s n (t), (1)
where v n and s n are members of some finite-dimensional spaces V h and S k composed of space and time functions having dimension N h and N k , respectively. We say that u m is an approximation of rank m of the exact solution, consisting of a summation of rank-one terms. We employ a greedy rank-one algorithm for computing the sequence of approximations. Specifically, the approximate solution u m is constructed iteratively, i.e., once u m-1 for m ≥ 1 is known (we set conventionally u 0 = 0), the functions v m and s m are computed by solving successively a global problem in space and in time having size N h and N k , respectively, and which is defined by minimizing some quadratic convex functional. The present method has the potential to be computationally effective if the exact solution can be approximated accurately by low-rank space-time tensors. In this case, space-time compression is meaningful and full parallelism in time can be exploited by the global time solves leading to an overall computational cost of the order of m(N h + N k ), whereas traditional time-stepping methods are expected to exhibit a cost of the order of N h × N k . Thus, computational benefits are expected whenever m min(N h , N k ). We notice that in the literature, approximate solutions of the form [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF] are typically obtained by a Proper Generalized Decomposition (PGD) [START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF][START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF][START_REF] Falcó | A proper generalized decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF][START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF], which is a greedy algorithm [START_REF] Temlyakov | Greedy approximation[END_REF]. In the context of parabolic problems, the PGD strategy has been first introduced within the LATIN method in [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation[END_REF]. Theoretical convergence results have been obtained in different contexts, see, e.g., [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF][START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF][START_REF] Cancès | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF][START_REF] Falcó | Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces[END_REF][START_REF] Cancès | Greedy algorithms for high-dimensional eigenvalue problems[END_REF]. We leave the question of adaptivity to future work, i.e., the discrete spaces V h and S k are fixed a priori in what follows. Possible applications of the present method can be envisaged in optimal control problems constrained by parabolic evolution equations (see, e.g., [START_REF] Gunzburger | Space-time adaptive wavelet methods for optimal control problems constrained by parabolic evolution equations[END_REF]) and in parabolic evolution equations with random input data (see, e.g., [START_REF] Hoang | Sparse tensor Galerkin discretization of parametric and random parabolic PDEs-analytic regularity and generalized polynomial chaos approximation[END_REF]); in both cases indeed, global space-time approaches are important. We also mention the dynamical low-rank integrators from, e.g., [START_REF] Koch | Dynamical low-rank approximation[END_REF][START_REF] Lubich | A projector-splitting integrator for dynamical low-rank approximation[END_REF][START_REF] Kieri | Discretized dynamical low-rank approximation in the presence of small singular values[END_REF] which can be used for parabolic evolution equations with the rank referring to the space variables.

Our starting point is the well-posed formulation of the parabolic evolution equation at hand in the setting of space-time Hilbert-Bochner spaces. More precisely, following Lions and Magenes [25, p. 234], the trial space is X = L 2 (I; V ) ∩ H 1 (I; V ) and the test space is Z = L 2 (I; V ) × L, where I is the (non-empty, bounded) time interval and the separable real Hilbert spaces (V, L, V ) form a Gelfand triple, i.e., V → L ≡ L → V with densely defined embeddings. The prototypical example is the heat equation for which

V = H 1 0 (Ω), L = L 2 (Ω), V = H -1 (Ω),
where Ω is a bounded, Lipschitz, open subset of R d , d ≥ 1. More generally, we consider time-dependent linear (possibly non-selfadjoint) operators A(t) : V → V that are bounded and coercive, pointwise in time (a.e.). One important assumption for the present method to remain computationally effective is that the space-time operator and source term in the parabolic evolution equation admit a separated representation in space and in time with a relatively low rank, see Eq. [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] below. This assumption is met in practice for a wide range of problems coming from the engineering and applied sciences. In several situations, it is also possible to devise such low-rank separated representations with good accuracy using the Empirical Interpolation Method (EIM) [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF]. The above well-posed space-time formulation allows one to view the parabolic evolution problem as a Minimal Residual (MinRes) formulation, where the exact solution is the unique minimizer over the trial space X of the (square of the) dual norm of the residual in Z . The present approximation method is formulated as a Galerkin method for the MinRes formulation. More precisely, we look for the minimizer over a finite-dimensional subspace X hk ⊂ X of the (square of the) dual norm of the residual measured with respect to a space semi-discrete subspace Z h ⊂ Z. Here, X hk is a linear space generated using space-time tensor-products of elements of the finite-dimensional subspaces

V h ⊂ V and S k ⊂ H 1 (I) considered above, whereas Z h = (V h ⊗ L 2 (I)) × V h , i.e., the time variable is not discretized in Z h .
Let us put our work in perspective with the literature. The subject of numerical methods to approximate parabolic evolution equations is extremely rich. One important class of methods are the traditional time-stepping schemes which approximate the solution on a succession of time sub-intervals by marching along the positive time direction, see, e.g., [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF] for an overview. In the context of parabolic evolution equations, implicit schemes are often preferred to circumvent the classical CFL restriction on explicit schemes, which is of the form

δ k δ 2 h (where δ k ∼ N -1 k is the time-step and δ h ∼ N -1/d h
is the space-step), but at the expense of having to solve a large linear system of equations at each time-step. The cost of having to solve these systems sequentially has motivated the devising of parareal methods [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF][START_REF] Gander | Analysis of the parareal timeparallel time-integration method[END_REF] based on iterative corrections at all time sub-intervals simultaneously using the global space-time discrete system associated with time-stepping methods. This global space-time discrete system has also been central in the devising of space-time domain decomposition methods based on waveform relaxation [START_REF] Janssen | Multigrid waveform relaxation of spatial finite element meshes: the continuous-time case[END_REF][START_REF] Gander | Overlapping Schwarz waveform relaxation for the heat equation in n dimensions[END_REF][START_REF] Giladi | Space-time domain decomposition for parabolic problems[END_REF]. In contrast to the above approaches which do not make a direct use of the well-posed functional setting in space-time Hilbert-Bochner spaces, a space-time adaptive wavelet method for parabolic evolution problems was proposed and analyzed in [START_REF] Schwab | Space-time adaptive wavelet methods for parabolic evolution problems[END_REF], involving a rather elaborate construction of the wavelet bases. Simpler hierarchical wavelet-type tensor bases on a space-time sparse grid were also considered in [START_REF] Griebel | A sparse grid space-time discretization scheme for parabolic problems[END_REF] within a heuristic space-time adaptive algorithm, but without offering guaranteed a-priori stability, uniformly with respect to the discretization parameters. We also mention the recent work [START_REF] Neumüller | Time-parallel iterative solvers for parabolic evolution equations[END_REF] where the above functional setting for parabolic evolution problems is used to devise preconditioned time-parallel iterative solvers. Furthermore, PGD approximations based on a discrete MinRes formulation measured in the space-time Euclidean norm of the components in a basis of X hk have been devised and evaluated numerically in [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF], obtaining promising results on various model parabolic evolution problems.

More recently, space-time Petrov-Galerkin discretizations of parabolic evolution equations were proposed and analyzed in the PhD Thesis [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF] and in the related papers [START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF]. Therein, the same MinRes formulation is considered as in the present work at the continuous level, and the approximate solution is typically sought in the same space-time discrete space. There are, however, several salient differences between [START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF] and the present work. First, in [START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF], the dual residual norm of the discrete solution is measured with respect to a fully discrete space-time test space, leading to a Petrov-Galerkin formulation, whereas we consider a space semi-discrete test space, leading to a standard Galerkin formulation. The difference is that a careful design of the test space is necessary in the Petrov-Galerkin setting. Precise results in this direction have been obtained in [START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF]. Let us for instance notice that, for a constant-in-time and selfadjoint differential operator in space (e.g., for the heat equation), the Crank-Nicolson scheme (obtained using continuous piecewise affine time-functions in the trial space and piecewise constant time-functions in the test space) is only conditionally stable, with a stability constant degenerating with the parabolic CFL condition, whereas unconditional stability is achieved by further refining the time mesh used to build the discrete test space as shown in [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF]Sec. 5.2.3]. In contrast with this, the present formulation automatically inherits uniform stability with respect to the time-step size. In particular, Lemma 3.1 below shows well-posedness and Lemma 3.3 a quasi-optimal error bound with a constant independent on the time discretization. Nonetheless, the condition number of the discrete matrices should behave as O(N -2 k ), so that, as usual, the precision can be limited in practice by round-off errors. The second difference lies in the way the discrete system of linear equations is solved iteratively: we use a greedy algorithm to build a sequence of approximate solutions having the form (1), whereas (a generalization of) the LSQR algorithm [START_REF] Paige | LSQR: an algorithm for sparse linear equations and sparse least squares[END_REF] is used in [START_REF] Andreev | Space-time discretization of the heat equation[END_REF]. Third, we report numerical results on a larger set of model parabolic evolution problems, including in particular non-selfadjoint operators of advection-diffusion-type at moderate Péclet numbers. Finally, let us mention, as already observed in [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF][START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF], that the inner product with which we equip the discrete test space Z h plays the role of a preconditioner in the discrete system of linear equations resulting from the discrete MinRes formulation. Equipping Z h with the natural norm leads to the appearance of the inverse of the stiffness matrix in space. Herein, we explore numerically the effect of relaxing the use of this preconditioner by simply equipping the space-part of Z h with the Euclidean norm of the components in a basis of V h . This corresponds to the approach considered in [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF].

In Section 2, we specify the functional setting for parabolic evolution equations and the MinRes formulation. This formulation is the basis for the discrete MinRes Galerkin formulation devised in Section 3, where one key idea is the use of a space semi-discrete test space to measure the dual norm of the residual. In Section 4, we present the greedy algorithm we consider to obtain a low-rank approximation of the discrete solution. In Section 5, we present two other discrete MinRes formulations, one using a fully discrete Petrov-Galerkin setting as in [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF][START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF] and one using the same space semi-discrete setting as in Section 3 but equipping the test space with the above-mentioned Euclidean norm. These additional formulations are introduced for the purpose of performing numerical comparisons. Numerical results on various test cases are discussed in Section 6. Finally, conclusions are drawn in Section 7.

Minimal Residual formulation of parabolic evolution equations

In this section, we present the MinRes formulation of parabolic equations, which is at the heart of the tensor approximation method we propose later on.

Parabolic equations

The functional setting for parabolic equations is well understood (see, e.g., the textbooks by Lions and Magenes [25, p. 234], Dautray and Lions [8, p. 513], and Wloka [38, p. 376]). Let

V → L ≡ L → V (2)
be a Gelfand triple where V and L are separable real Hilbert spaces respectively equipped with inner products •, • V and •, • L , with associated norms • V and • L . The symbol → represents a densely defined and continuous embedding. Let T > 0 be the time horizon and let I := (0, T ) be the time interval. Let A : I → L(V, V ) be a strongly measurable time-function with values in the Hilbert space of bounded linear operators from V to V . We assume that the following boundedness and coercivity properties hold true: there exist 0 < α ≤ M < +∞ such that a.e. t ∈ I,

A(t)v V ≤ M v V , ∀v ∈ V, (3a)
A(t)v, v V ,V ≥ α v 2 V , ∀v ∈ V. (3b)
We do not require A(t) to be selfadjoint. Let us define the Hilbert-Bochner spaces

X := L 2 (I; V ) ∩ H 1 (I; V ), Y := L 2 (I; V ), Z := Y × L. (4)
Since X → C 0 (I; L) with I = [0, T ], the value at any time t ∈ I of any function x ∈ X is well-defined as an element of L. In particular, we denote x(0) ∈ L the initial value of x at the time t = 0. The spaces X and Z are equipped with the norms

x 2 X := x 2 L 2 (I;V ) + M -2 ∂ t x 2 L 2 (I;V ) + α -1 x(T ) 2 L , ∀x ∈ X, (5a)
z 2 Z := y 2 L 2 (I;V ) + α -1 g 2 L , ∀z = (y, g) ∈ Y × L = Z, (5b)
where the various scaling factors are introduced to be dimensionally consistent. Let f ∈ Y = L 2 (I; V ) and let u 0 ∈ L. We consider the following parabolic problem: find u ∈ X such that

∂ t u(t) + A(t)u(t) = f (t), in V , a.e. t ∈ I, u(0) = u 0 , in L. (6)
For the present space-time tensor methods to be computationally effective, we assume that the operator A and the source term f have the following separated form:

A(t) = 1≤p≤P µ (p) (t)A (p) ∈ L(V, V ), f (t) = 1≤q≤Q λ (q) (t)f (q) ∈ V , (7)
for some positive integers P, Q taking moderate values, where µ (p) ∈ L ∞ (I), A (p) ∈ L(V, V ) for all 1 ≤ p ≤ P , and λ (q) ∈ L 2 (I), f (q) ∈ V for all 1 ≤ q ≤ Q. A similar decomposition is considered, e.g., in [START_REF] Mantzaflaris | Low-rank space-time isogeometric analysis for parabolic problems with varying coefficients[END_REF] for the space-time isogeometric discretization of parabolic problems with varying coefficients.

Example 2.1 (Heat equation).

Let Ω be a Lipschitz domain in R d , V = H 1 0 (Ω), L = L 2 (Ω), and V = H -1 (Ω). Let µ : I → R be a measurable function bounded from above and from below away from zero uniformly on I. Then, the time-dependent family of operators (A(t)) t∈I defined such that A(t) := -µ(t)∆, for a.e. t ∈ I, where ∆ ∈ L(V ; V ) is the Laplacian operator, satisfies the assumptions (3) and [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF]. Whenever µ(t) ≡ 1, the family of operators is timeindependent, and one recovers the prototypical example of the heat equation.

Well-posedness and minimal residual formulation

It is convenient to introduce the operator A : X → Y so that (Ax)(t) = ∂ t x(t) + A(t)x(t) ∈ V , a.e. t ∈ I. (8)
Problem ( 6) can be equivalently rewritten using the operator B :

X → Z = Y × L such that Bx = (Ax, x(0)) ∈ Z , ∀x ∈ X. (9)
Then, an equivalent reformulation of problem (6) reads as follows:

find u ∈ X such that Bu = (f, u 0 ) in Z . (10)
It is well-known that the operator B is bounded, i.e., B ∈ L(X, Z ), and satisfies the following two properties:

∃β > 0 s.t. inf x∈X sup z∈Z Bx, z Z ,Z x X z Z ≥ β, (11a)
∀z ∈ Z, ( Bx, z Z ,Z = 0, ∀x ∈ X) =⇒ (z = 0), ( 11b
)
where it is implicitly understood that nonzero arguments are considered in the inf-sup condition, and where, for all (y, g)

∈ Y × L = Z, Bx, (y, g) Z ,Z = Ax, y Y ,Y + x(0), g L = ˆI ∂ t x(t) + A(t)x(t), y(t) V ,V dt + x(0), g L . (12)
Therefore, owing to the Banach-Nečas-Babuška Theorem (see, e.g., [9, Thm. 2.6]), B is an isomorphism. The proof of the well-posedness of parabolic problems by means of inf-sup arguments can be found in [9, Thm. 6.6] using a strongly enforced initial condition; a systematic treatment can be found more recently in [START_REF] Tantardini | The L 2 -projection and quasi-optimality of Galerkin methods for parabolic equations[END_REF]. Since the operator norm of B and the inf-sup constant β in (11a) play an important role in what follows, we provide respectively an upper bound and a lower bound for these two constants. For a.e. t ∈ I, the inverse adjoint operator A(t) -T is well-defined in L(V ; V ), and we have

M -1 ϕ V ≤ A(t) -T ϕ V ≤ α -1 ϕ V and ϕ, A(t) -T ϕ V ,V ≥ α M 2 ϕ 2 V , for all ϕ ∈ V . Lemma 2.2 (Boundedness). The norm of the operator B : X → Z is such that B L(X;Z ) := sup x∈X sup z∈Z Bx, z Z ,Z x X z Z ≤ √ 3M. Proof. For all x ∈ X, z = (y, g) ∈ Z = Y × L, we have Bx, z Z ,Z = Ax, y Y ,Y + x(0), g L ≤ Ax Y y Y + x(0) L g L ≤ Ax 2 Y + α x(0) 2 L 1/2 y 2 Y + α -1 g 2 L 1/2 ≤ √ 3M x X z Z , since Ax 2 Y ≤ 2M 2 (M -2 ∂ t x 2 Y + x 2 Y ) ≤ 2M 2 x 2 X and α x(0) 2 L ≤ α x(T ) 2 L +2α ∂ t x Y x Y ≤ α x(T ) 2 L +αM (M -2 ∂ t x 2 Y + x 2 Y ) ≤ M 2 x 2 X (recall that α ≤ M ). Lemma 2.3 (Inf-sup constant). (11a) holds true with the inf-sup constant β ≥ α 2 M (1 + κ 2 ) -1 2
, where κ := ess sup t∈I

1 2 A(t)A(t) -T -I L(V ;V ) . Proof.
For completeness, we briefly outline the proof which is similar to that of [START_REF] Tantardini | The L 2 -projection and quasi-optimality of Galerkin methods for parabolic equations[END_REF]Prop. 3.1], but with a different scaling for the norms. Let x ∈ X and let us take

z = (A(t) -T ∂ t x(t) + x(t), x(0)) ∈ Z.
Then, using the coercivity of A(t) and of A(t) -1 , we have

Bx, z Z ,Z = ˆI ∂ t x(t) + A(t)x(t), A(t) -T ∂ t x(t) + x(t) V ,V dt + x(0) 2 L ≥ αM -2 ∂ t x 2 L 2 (I;V ) + α x 2 L 2 (I;V ) + x(T ) 2 L ≥ α x 2 X .
Moreover, using again the coercivity of A(t) and the boundedness of A(t) and A(t) -T , we have

z 2 Z = ˆI A(t) -T ∂ t x(t) + x(t) 2 V dt + α -1 x(0) 2 L ≤ α -1 ˆI A(t) A(t) -T ∂ t x(t) + x(t) , A(t) -T ∂ t x(t) + x(t) V ,V dt + α -1 x(0) 2 L ≤ M α -1 x 2 L 2 (I;V ) + M 2 α -2 ∂ t x 2 L 2 (I;V ) + α -1 x(T ) 2 L + 2κα -1 ∂ t x L 2 (I;V ) x L 2 (I;V ) ≤ (1 + κ 2 )M 2 α -2 x 2 X
, and the conclusion is straightforward.

Remark 2.4 (Heat equation). Sharp estimates of the inf-sup constant β for the heat equation (with µ(t) ≡ 1 on I so that α = M = 1 and κ = 0) can be found in [START_REF] Urban | A new error bound for reduced basis approximation of parabolic partial differential equations[END_REF][START_REF] Ern | Locally Space-Time Efficient, and Polynomial-Degree Robust a Posteriori Error Estimates for High-Order Discretizations of Parabolic Problems[END_REF] using the above norms.

The solution to the parabolic equation ( 6) is the global minimizer of the residual-based quadratic functional E : X → R, defined such that

E(x) := 1 2 Bx -(f, u 0 ) 2 Z , ∀x ∈ X, (13)
where we equip the space

Z = Y × L with the norm (φ, g) 2 Z := φ 2 L 2 (I;V ) + α g 2 L . (14)
Since the functional E is strongly convex on X with parameter β 2 > 0 owing to the inf-sup condition (11a), E admits a unique global minimizer in X, and since the operator B is surjective, the minimum value of E on X is zero. In other words, the unique solution to (6) can be equivalently characterized as follows:

u = argmin x∈X 1 2 Ax -f 2 Y + α x(0) -u 0 2 L . ( 15
)
3 Discrete Minimal Residual Galerkin formulation

In this section, we introduce a discrete energy based on a space semi-discrete Galerkin method to approximate the unique minimizer in [START_REF] Giladi | Space-time domain decomposition for parabolic problems[END_REF]. We consider finite-dimensional spaces V h and S k such that

V h ⊂ V, S k ⊂ S := H 1 (I), (16)
and we set N h := dim(V h ) and N k := dim(S k ). Typically, V h is constructed using P 1 Lagrange finite elements on a space mesh of Ω and S k is constructed using continuous, piecewise affine functions on a time mesh of I. We are going to seek the discrete minimizer in the tensor-product space

X hk := V h ⊗ S k ⊂ X, (17)
which is of dimension dim(X hk ) = N h × N k .

Space semi-discrete Galerkin approximation

Let us set

X h = V h ⊗H 1 (I), Y h := V h ⊗L 2 (I) ≡ L 2 (I; V h ), Y h = V h ⊗L 2 (I) ≡ L 2 (I; V h ). (18) Since V h is a subspace of V and owing to (2), we have Y h ⊂ Y and Y ⊂ Y h ,
where the second inclusion follows by restricting the action of linear forms on

V to V h . Let us define f h ∈ Y h s.t. f h (t) = f (t)| V h a.e. t ∈ I.
Recalling the separated form [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF], we have

f h (t) := 1≤q≤Q λ (q) (t)f (q) | V h a.e. t ∈ I. ( 19
) Let u 0h be the L-orthogonal projection of u 0 onto V h . Let A h : X h → Y h be s.t. (A h x h )(t) := (J V h ⊗ ∂ t )x h (t) + A h (t)x h (t) ∈ V h , a.e. t ∈ I, (20)
where

J V h : V h → V h is the injection resulting from (2), i.e., J V h v h , w h V h ,V h = v h , w h L for all v h , w h ∈ V h , and A h (t) : V h → V h is the discrete counterpart of A(t) s.t. A h (t)v h , w h V h ,V h = A(t)v h , w h V ,V . Let us introduce the space semi-discrete space Z h = Y h × L h so that Z h = Y h × L h , where L h coincides with V h as linear space but is equipped with the norm of L (note that L h ⊂ L ≡ L ⊂ L h ). Note that Z h ⊂ Z. Let B h : X h → Z h be the operator defined for x h ∈ X h by B h x h = (A h x h , x h (0)) ∈ Z h = Y h × L h , (21)
and such that, for all (y h , g h ) ∈ Y h × L h , (compare with (12))

B h x h , (y h , g h ) Z h ,Z h = A h x h , y h Y h ,Y h + x h (0), g h L = Bx h , (y h , g h ) Z ,Z . (22)
The space semi-discrete formulation is as follows: find u h ∈ X h such that

B h u h = (f h , u 0h ) in Z h . (23)
The well-posedness of this formulation is ensured by the following lemma, where the subspace X h ⊂ X is equipped with the norm x h 2

X h := x h 2 L 2 (I;V ) + M -2 ∂ t x h 2 L 2 (I;V h ) + α -1 x h (T ) 2
L , and the subspace

Z h ⊂ Z is equipped with the norm of Z. Recall the inf-sup constant β = α 2 M (1 + κ 2 ) -1 2 from Lemma 2.3. Lemma 3.1. The operator B h : X h → Z h satisfies B h L(X h ;Z h ) ≤ √ 3M, and
inf x h ∈X h sup z h ∈Z h B h x h , z h Z h ,Z h x h X h z h Z ≥ β, ( 24a
)
∀z h ∈ Z h , ( B h x h , z h Z h ,Z h = 0, ∀x h ∈ X h ) =⇒ (z h = 0). ( 24b
)
Proof. The upper bound on B h L(X h ;Z h ) is shown as in the proof of Lemma 2.2.

To prove (24a), one can use the same arguments as in the proof of Lemma 2.3 by picking in the supremizing set

z h = (A h (t) -T ∂ t x h (t) + x h (t), x h (0)) ∈ Y h × L h = Z h .
Finally, one can prove (24b) using the following arguments, as in [9, Thm. 6.6]. Let

z h = (y h , g h ) ∈ Z h . Taking x h arbitrary in V h ⊗ C ∞ 0 (I) shows that (J V h ⊗ ∂ t )y h -A -T h y h = 0 in Y h .
Taking next x h = tv h with v h arbitrary in V h proves that y h (T ) = 0, and taking x h = ty h , one concludes that y h = 0. Finally, taking

x h = v h with v h arbitrary in V h yields g h = 0. Lemma 3.1 implies that B h : X h → Z h is an isomorphism such that β x h X h ≤ B h x h Z h ≤ √ 3M x h X h , (25)
for all x h ∈ X h . The solution u h to the equation ( 23) is the unique minimizer of the discrete energy functional E h : X h → R defined for all x h ∈ X h by

E h (x h ) := 1 2 B h x h -(f h , u 0h ) 2 Z h = 1 2 A h x h -f h 2 Y h + α x h (0) -u 0h 2 L . (26
) An important property of this discrete energy functional is the strong convexity that is inherited from the continuous setting, uniformly with respect to the space discretization parameter. More precisely, the functional E h is strongly convex on X h with parameter β 2 > 0 owing to the inf-sup condition (24a). Since the operator B h is surjective, the minimum value of E h on X h is zero and is attained at u h .

Remark 3.2 (Norm • X h ). The difference between the • X -norm and the

• X h -norm lies in the use of the dual norm in V h and not in V to measure the time-derivative. Note that x h X h ≤ x h X , for all x h ∈ X h . The reason for this difference is that, as shown in [START_REF] Tantardini | The L 2 -projection and quasi-optimality of Galerkin methods for parabolic equations[END_REF], the equivalence of the two norms, uniformly with respect to the space discretization, holds true if and only if the L-orthogonal projection onto V h is V -stable. This uniform stability (with V = H 1 0 (Ω) and L = L 2 (Ω)) is, in turn, not known to hold true if general shaperegular meshes are used to build the finite element space V h ; it does hold true if quasi-uniform meshes are used (as it is the case in the present numerical experiments). We emphasize that the use of a discrete dual norm to measure the time-derivative is a general feature that arises in the quasi-optimality of space semi-discrete Galerkin methods for parabolic evolution problems [START_REF] Tantardini | The L 2 -projection and quasi-optimality of Galerkin methods for parabolic equations[END_REF], and is not specific to the present setting.

Minimal residual Galerkin approximation

An approximation u hk ∈ X hk of u h ∈ X h is now defined as the unique minimizer of the discrete energy functional E h restricted to the subspace X hk of X h , i.e.we look for

u hk = argmin x hk ∈X hk E h (x hk ) = argmin x hk ∈X hk 1 2 A h x hk -f h 2 Y h + α x hk (0) -u 0h 2 L . ( 27
)
We emphasize the use of the space semi-discrete test space Y h to measure the dual norm of the residual. We obtain the following quasi-optimal error estimate. Lemma 3.3 (Error estimate). Let u h be the unique solution to [START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF], and let u hk ∈ X hk be the unique minimizer of [START_REF] Mantzaflaris | Low-rank space-time isogeometric analysis for parabolic problems with varying coefficients[END_REF]. Then, we have

u h -u hk X h ≤ C inf x hk ∈X hk u h -x hk X h , (28)
where

C = √ 3M β
is independent of the time discretization.

Proof. Using (25), we have

u hk -u h X h ≤ β -1 B h (u hk -u h ) Z h = β -1 min x hk ∈X hk B h (x hk -u h ) Z h ≤ √ 3M β -1 min x hk ∈X hk x hk -u h X h ,
which proves the assertion.

The unique minimizer of the quadratic discrete minimization problem (27) can be characterized by a system of linear equations. To write this system, let us first introduce the Riesz isomorphism

R V h : V h → V h such that R V h v h , w h V h ,V h = v h , w h V for all v h , w h ∈ V h (note that R V h differs from the injection J V h in- troduced above). Let R Y h : Y h → Y h be the space-time Riesz isomorphism such that R Y h = R V h ⊗ I L 2 , ( 29
)
where I L 2 is the identity operator in L 2 (I) (it is actually the Riesz isomorphism from L 2 (I) onto L 2 (I) ≡ L 2 (I)). The quadratic discrete minimization problem ( 27) is equivalent to the following linear problem: find u hk ∈ X hk such that

B hk u hk = g hk , (30)
with B hk : X hk → X hk and g hk ∈ X hk such that, for all x hk , z hk ∈ X hk ,

B hk x hk , z hk X hk ,X hk = A h x hk , R -1 Y h A h z hk Y h ,Y h + α x hk (0), z hk (0) L , (31a) g hk , z hk X hk ,X hk = f h , R -1 Y h A h z hk Y h ,Y h + α u 0h , z hk (0) L . ( 31b
)
Let us briefly describe the algebraic realization of the discrete problem [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF]. Let (ψ i ) 1≤i≤N h be a basis of V h and let (φ l ) 1≤l≤N k be a basis of S k . We can then seek for the components of the unique solution u hk of (30) in the basis

(ψ i ⊗ φ l ) 1≤i≤N h ,1≤l≤N k of X hk , i.e., we seek u = (u il ) 1≤i≤N h ,1≤l≤N k ∈ R N h N k such that u hk = N h i=1 N k l=1 u il ψ i ⊗ φ l . (32)
We define the following matrices of size N h ×N h (related to the space discretization):

(D h ) ij = ψ j , ψ i V , (M h ) ij = ψ j , ψ i L , (33)
and the following matrices of size N k × N k (related to the time discretization):

(D k ) lm = ˆI φ m (t)φ l (t)dt, (M k ) lm = ˆI φ m (t)φ l (t)dt, (O k ) lm = φ m (0)φ l (0). ( 34
)
Recalling the separated form [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF], we introduce the following matrix of size

N h × N h : (A (p) h ) ij = A (p) ψ j , ψ i V ,V , (35)
and the following matrices of size

N k × N k : (M (p,p ) k ) lm = ˆI µ (p) (t)µ (p ) (t)φ m (t)φ l (t)dt, (E (p) k ) lm = ˆI µ (p) (t)φ m (t)φ l (t)dt, (36)
for all 1 ≤ p, p ≤ P . Then, we obtain the following symmetric positive-definite linear system in

R N h N k : Bu = g, (37)
with the matrix

B = M h D -1 h M h ⊗ D k + 1≤p≤P 2sym (A (p) h ) T D -1 h M h ⊗ E (p) k + 1≤p,p ≤P (A (p) h ) T D -1 h A (p ) h ⊗ M (p,p ) k + αM h ⊗ O k , (38)
where sym(

Z h ⊗ Z k ) = 1 2 (Z h ⊗ Z k + Z T h ⊗ Z T k )
for any matrix Z h of size N h × N h and any matrix Z k of size N k × N k , and the right-hand side

g = 1≤q≤Q M h D -1 h f (q) h ⊗ e (q) k + 1≤p≤P 1≤q≤Q (A (p) h ) T D -1 h f (q) h ⊗ d (p,q) k + αu 0h ⊗ i k , (39)
with the vectors (f

(q) h ) i = f (q) | V h , ψ i V ,V = f (q) , ψ i V ,V , (u 0h ) i = u 0h , ψ i L = u 0 , ψ i L , for all 1 ≤ i ≤ N h , and (e (q) k ) l = ´I λ (q) (t)φ l (t)dt, (d (p,q) k ) l = ´I µ (p) (t)λ (q) (t)φ l (t)dt, (i k ) l = φ l (0), for all 1 ≤ l ≤ N k .
Example 3.4 (Heat equation). Let us consider the heat equation where P = 1, µ (1) (t) ≡ 1 and A (1) = -∆, and let us equip the space V with the H 1 -seminorm so that v, w V = ´Ω ∇v(x)•∇w(x)dx = A (1) v, w V ,V . Then the expression of B simplifies as follows:

B = M h D -1 h M h ⊗ D k + M h ⊗ 2sym(E k ) + D h ⊗ M k + αM h ⊗ O k , ( 40
)
with the following matrix of size N k × N k :

(E k ) lm = ˆI φ m (t)φ l (t)dt. ( 41
)
4 Low-rank approximation

In this section, we present the low-rank approximation method we use to approximate iteratively the unique minimizer of ( 27) (or, equivalently, the unique solution to the linear system (37)). We consider here a greedy algorithm [START_REF] Temlyakov | Greedy approximation[END_REF][START_REF] Cancès | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF][START_REF] Falcó | Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces[END_REF][START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] which is an iterative procedure such that, at each iteration m ∈ N * , one computes an approximation u m hk ∈ X hk of the solution u hk ∈ X hk of ( 27) in the form

u m hk (x, t) = 1≤n≤m v n h (x) ⊗ s n k (t), ( 42
) with v n h ∈ V h and s n k ∈ S k for all 1 ≤ n ≤ m.
The algorithm can be outlined as follows:

GREEDY ALGORITHM:

1. Set u 0 hk = 0 and m = 1.

Solve for

(v m h , s m k ) ∈ argmin (v h ,s k )∈V h ×S k E h (u m-1 hk + v h ⊗ s k ). ( 43
) Set u m hk := u m-1 hk + v m h ⊗ s m k .
3. Check convergence, and if not satisfied, set m ← m + 1 and go to step (2).

The following relative stagnation-based stopping criterion is used with a tolerance greedy > 0:

v m h ⊗ s m k X u m hk X < greedy . (44)
Using the general results from [START_REF] Cancès | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF][START_REF] Falcó | Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces[END_REF], one can verify that the iterations of the above greedy algorithm are well-defined using the discrete minimal residual formulation presented in Section 3. Recall that the uniqueness of the solution to the minimization problem [START_REF] Mantzaflaris | Low-rank space-time isogeometric analysis for parabolic problems with varying coefficients[END_REF] follows from the strong convexity of the functional E h , and the sequence (u m hk ) m∈N converges to u hk as n goes to infinity. Actually, it can be checked that this convergence result still holds true in the infinite-dimensional setting.

In the above greedy algorithm, the minimization problem (43) is nonlinear. Therefore, it is not straightforward to solve it and in practice, one often considers an alternating minimization algorithm (see [START_REF] Uschmajew | Local convergence of the alternating least squares algorithm for canonical tensor approximation[END_REF]), based on the following fixedpoint iterative scheme: ALTERNATING MINIMIZATION ALGORITHM FOR (43):

1. Choose s m,0 k ∈ S k randomly and set p = 1. 2. Let v m,p h ∈ V h be the unique solution to v m,p h = argmin v h ∈V h E h u m-1 hk + v h ⊗ s m,p-1 k . ( 45
)
Compute s m,p k ∈ S k to be the unique solution to

s m,p k = argmin s k ∈S k E h u m-1 hk + v m,p h ⊗ s k . (46)
3. Check convergence, and if not satisfied, set p ← p + 1 and go to step (2).

The following relative stagnation-based stopping criterion is used with a tolerance alt > 0:

v m,p h ⊗ s m,p k -v m,p-1 h ⊗ s m,p-1 k X v m,p h ⊗ s m,p k X < alt . ( 47
)
The cost of an iteration of the alternating minimization algorithm is of order (N h + N k ). Provided the number of fixed-point iterations remains reasonably small, the cost of each iteration of the greedy algorithm can be estimated to scale also as (N h + N k ). We will verify in our numerical experiments that this is indeed the case.

(v m h , s m k ) = argmin (v h ,s k )∈R N h ×R N k 1 2 u m-1 hk + v h ⊗ s k T B u m-1 hk + v h ⊗ s k -u m-1 hk + v h ⊗ s k T g ,
where u m-1 hk denotes the vector in R N h N k containing the coordinates of u m-1 hk in the basis (ψ i ⊗ φ l ) 1≤i≤N h ,1≤l≤N k . Similarly, the matrix form of problems (45) and (46) is as follows:

v m,p h = argmin v h ∈R N h 1 2 u m-1 hk + v h ⊗ s m,p-1 k T B u m-1 hk + v h ⊗ s m,p-1 k -u m-1 hk + v h ⊗ s m,p-1 k T g , s m,p k = argmin s k ∈R N k 1 2 u m-1 hk + v m,p h ⊗ s k T B u m-1 hk + v m,p h ⊗ s k -u m-1 hk + v m,p h ⊗ s k T g .

Other discrete minimal residual methods

In this section, we describe for the purpose of numerical comparisons in Section 6 two other discrete minimal residual approaches. The discrete method introduced in Section 3 is henceforth referred to as Method 1. The first variant, called Method 2, hinges on a fully discrete Petrov-Galerkin setting as devised in [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF][START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF]. The second variant, called Method 3, uses the same space semi-discrete setting as Method 1, but the test space is now equipped with a simple Euclidean norm of the components on a basis of V h ; Method 3 has been introduced in [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF].

Method 2: fully discrete Petrov-Galerkin method

Let us set

Y hk := V h ⊗ S P k , Y hk = V h ⊗ (S P k ) , (48)
where

S P k is a finite-dimensional subspace of L 2 (I) so that S P k ⊂ L 2 (I) ≡ L 2 (I) ⊂ (S P k ) . The injection J S P k : L 2 (I) → (S P k ) is such that J S P k = R (S P k ) • Π S P k where Π S P k is the L 2 (I)-orthogonal projection from L 2 (I) onto S P k and R (S P k ) : S P k → (S P k )
is the Riesz isomorphism so that R (S P k ) q, r (S P k ) ,S P k = q, r L 2 (I) , for all q, r ∈ S P k . Let us set dim(S P k ) = N P k . Recalling the separated form [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF], let us define f hk ∈ Y hk s.t.

f hk (t) = 1≤q≤Q (J S P k λ (q) )(t)f (q) | V h a.e. t ∈ I. (49)
Let A hk := (I V h ⊗ J S P k )A h : X hk → Y hk where I V h is the identity operator in V h and A h is defined by [START_REF] Kieri | Discretized dynamical low-rank approximation in the presence of small singular values[END_REF]. We consider the discrete energy functional E fdPG hk : X hk → R defined as

E fdPG hk (x hk ) := 1 2 A hk x hk -f hk 2 Y hk + α x hk (0) -u 0h 2 L , ∀x hk ∈ X hk .
(50) The discrete minimization problem is as follows: find u fdPG hk ∈ X hk such that

u fdPG hk = argmin x hk ∈X hk E fdPG hk (x hk ). ( 51
)
Remark 5.1 (Comparison of energies). Since f hk = (I V h ⊗ J S P k )f h with f h defined by [START_REF] Janssen | Multigrid waveform relaxation of spatial finite element meshes: the continuous-time case[END_REF], we have

A hk x hk -f hk Y hk = (I V h ⊗ J S P k )(A h x hk -f h ) Y hk ≤ A h x hk -f h Y h , which implies that E fdPG hk (x hk ) ≤ E h (x hk ) for all x hk ∈ X hk .
As shown in [START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF] in the case of time-independent and selfadjoint operators A ∈ L(V ; V ), the Hessian of the discrete energy E fdPG hk induces a bilinear form that satisfies an inf-sup condition that degenerates with the parabolic CFL. In the general case with a time-dependent differential operator, the positivity of the inf-sup constant is not guaranteed a priori, which means that the discrete energy functional E fdPG hk may be only convex in some unfavorable situations. This means that in such cases, global minimizers of (51) exist but may not be unique. Any minimizer satisfies the following system of linear equations: find

u fdPG hk ∈ X hk such that B fdPG hk u fdPG hk = g fdPG hk , (52)
with B fdPG hk : X hk → X hk and g fdPG hk ∈ X hk such that, for all x hk , z hk ∈ X hk ,

B fdPG hk x hk , z hk X hk ,X hk = A hk x hk , R -1 Y hk A hk z hk Y hk ,Y hk + α x hk (0), z hk (0) L , (53a)
g fdPG hk , z hk X hk ,X hk = f hk , R -1 Y hk A hk z hk Y hk ,Y hk + α u 0h , z hk (0) L , (53b)
with the space-time Riesz isomorphism

R Y hk = R V h ⊗ J S P k : Y hk → Y hk . Let us point out that, since J S P k = R (S P k ) on S P k , we have R Y hk = R V h ⊗ R (S P k )
. Let us briefly describe the algebraic realization of the discrete problem (52). Recall that (ψ i ) 1≤i≤N h is a basis of V h and (φ l ) 1≤l≤N k is a basis of S k . Let (φ P l ) 1≤l≤N P k be a basis of S P k . In addition to the square matrices M h , D h , O k defined by ( 33), ( 34), [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF], we consider the square matrix M P k of size N P k × N P k and the rectangular matrix

E PG k of size N P k × N k such that (M P k ) lm = ˆI φ P m (t)φ P l (t)dt, (E PG k ) lm = ˆI φ m (t)φ P l (t)dt, ( 54
)
and the rectangular matrices M

PG,(p) k

, for all 1 ≤ p ≤ P , of size

N P k × N k such that (M PG,(p) k ) lm = ˆI µ (p) (t)φ m (t)φ P l (t)dt. ( 55
)
Then, we obtain the following symmetric positive-definite linear system in R N h N k :

B fdPG u fdPG = g fdPG , (56)
with the matrix

B fdPG = M h D -1 h M h ⊗ (E PG k ) T (M P k ) -1 E PG k + 1≤p≤P 2sym (A (p) h ) T D -1 h M h ⊗ (M PG,(p) k ) T (M P k ) -1 E PG k + 1≤p,p ≤P (A (p) h ) T D -1 h A (p ) h ⊗ (M PG,(p) k ) T (M P k ) -1 M PG,(p ) k + αM h ⊗ O k , (57)
and the right-hand side

g fdPG = 1≤q≤Q M h D -1 h f (q) h ⊗ (E PG k ) T (M P k ) -1 e P,(q) k + 1≤p≤P 1≤q≤Q (A (p) h ) T D -1 h f (q) h ⊗ (M PG,(p) k ) T (M P k ) -1 d P,(q) k + αu 0h ⊗ i k , (58)
with (e P,(q) k ) l = ´I λ (q) (t)(φ P l ) (t)dt and (d P,(q) k ) l = ´I λ (q) (t)φ P l (t)dt, for all 1 ≤ l ≤ N P k . Remark 5.2 (Lowest-order Petrov-Galerkin discretization). Assume that S k is composed of continuous, piecewise affine functions and that S P k is composed of piecewise constant functions on the same time mesh so that dim(S

P k ) = N k -1.
This corresponds to the well-known Crank-Nicolson time scheme. Then, one can readily verify that

(E PG k ) T (M P k ) -1 E PG k = D k , (M PG,(p) k ) T (M P k ) -1 E PG k = E (p) k , (59)
for all 1 ≤ p ≤ P , where D k is defined in [START_REF] Temlyakov | Greedy approximation[END_REF] and [START_REF] Urban | A new error bound for reduced basis approximation of parabolic partial differential equations[END_REF]. As a consequence, there is only one term composing the matrices B and B fdPG that differs, namely the time matrix in the double summation over p, p where this matrix is

E (p) k in (
M (p,p ) k for B and is (M PG,(p) k ) T (M P k ) -1 M PG,(p ) k
for B fdPG . Even for the heat equation with P = 1 and µ (1) (t) ≡ 1, these matrices, which become, respectively,

M k and (M PG k ) T (M P k ) -1 M PG k with (M PG k ) lm = ´I φ m (t)φ P l (t)dt, are still dif- ferent. Note that M k ≥ (M PG k ) T (M P k ) -1 M PG k
in the sense of quadratic forms, which is compatible with our above observation on the discrete energies that E fdPG hk (x hk ) ≤ E h (x hk ) for all x hk ∈ X hk . As observed in [START_REF] Andreev | Stability of space-time Petrov-Galerkin discretizations for parabolic evolution equations[END_REF][START_REF] Andreev | Stability of sparse space-time finite element discretizations of linear parabolic evolution equations[END_REF][START_REF] Andreev | Space-time discretization of the heat equation[END_REF], uniform stability with respect to the time discretization is not guaranteed, but this can be fixed, e.g., by constructing the discrete test space S P k using a time mesh that is twice as fine as that used for the discrete trial space.

Method 3: an unpreconditioned space semi-discrete Galerkin method

We consider the same space semi-discrete setting as in Section 3 but we now equip the space V h with the Euclidean norm of the components on the basis (ψ i ) 1≤i≤N h instead of considering as before the norm induced by V . The main motivation for this change is that it avoids the appearance of the inverse stiffness matrix D -1 h in the linear system. We obtain the following symmetric positivedefinite linear system in R N h N k :

B unprec u unprec = g unprec , ( 60
)
with the matrix

B unprec = (M h I h M h ) ⊗ D k + 1≤p≤P 2sym ((A (p) h ) T I h M h ) ⊗ E (p) k + 1≤p,p ≤P ((A (p) h ) T I h A (p ) h ) ⊗ M (p,p ) k + αM h ⊗ O k , (61)
and the right-hand side

g unprec = 1≤q≤Q M h I h f (q) h ⊗e (q) k + 1≤p≤P 1≤q≤Q (A (p) h ) T I h f (q) h ⊗d (p,q) k +αu 0h ⊗i k , ( 62
)
where I h is the identity matrix of size N h × N h . The present formulation is chosen for illustrative purposes; in practice, one can also replace D -1 h by another matrix.

Numerical results

For all the test cases, we consider the space domain Ω = (0, 1) × (0, 1), the time interval I = [0, 1], and the functional spaces V = H 1 0 (Ω) and L = L 2 (Ω).

We consider first the heat equation, where the differential operator A is timeindependent and selfadjoint, then we consider a time-oscillatory diffusion problem, where the operator is time-dependent and selfadjoint, and finally a convectiondiffusion equation, where the operator is time-independent and non-selfadjoint. The scaling factor for the contribution of the initial condition to the residual functional is always taken to be α := 1. Let T h be a shape-regular mesh of the domain Ω; in what follows, we consider uniform meshes composed of square cells. The finite-dimensional finite element subspace V h ⊂ V of dimension N h is composed of continuous, piecewise bilinear functions on T h vanishing at the boundary. Let T k be a mesh of the interval I; for simplicity, we consider uniform meshes in time. The finite-dimensional subspace S k ⊂ H 1 (I) of dimension N k is composed of continuous, piecewise affine functions on T k . In what follows, all the norms of residuals of algebraic quantities are evaluated using the Euclidean norm in R N h N k which is denoted • 2 . When comparing to Method 2 (see Section 5.1), we considered the Crank-Nicolson time scheme discussed in Remark 5.2. We also performed systematic comparisons with the uniformly stable variant using a finer time mesh for the test space, but we did not observe any significant difference in the results obtained for all the test cases considered herein.

Test case 1: heat equation with manufactured solution

We consider the heat equation with the time-independent, selfadjoint operator A = -∆. The initial condition is zero and the source term is evaluated from the following manufactured solution:

u(x, y, t) = 1≤n≤10 n -4 sin(πn 3 t)sin(πnx)sin(πny).

The discretization parameters are N h = (2 6 ) 2 and N k = 2 13 , and the stopping tolerances are greedy = 10 -5 and alt = 5 × 10 -2 . The left panel of Figure 1 presents the decrease of the relative residual as a function of the number of greedy iterations for Methods 1, 2 and 3. More precisely, we plot B i u m i -g i 2 / g i 2 where i ∈ {1, 2, 3} is the method index and m is the greedy iteration counter. We notice that the three methods take about the same number of iterations (13, 15, and 17, respectively). This number is slightly larger than the space-time rank of the manufactured exact solution which is equal to 10. However, the relative residual takes larger values for Method 3 than for Methods 1 and 2. The right panel of Figure 1 presents the cumulated number of alternating minimization iterations in the greedy algorithm for Methods 1, 2 and 3. We observe that this number is about the same for Methods 1 and 2, whereas it is about 1.8 times larger for Method 3. Therefore, the use of a preconditioner, although it requires some additional computational effort, is beneficial to the efficiency of the overall behavior of the greedy algorithm. It is interesting to notice that with Methods 1 and 2, we have solved at convergence of the greedy algorithm about 50 linear systems in space, which is about 0.25% of the amount that would have been solved by using an implicit time-stepping method (recall that N k = 2 13 ). We can make two further remarks concerning the decrease of the (relative) residual. First, we can decompose the residual of the space-time linear system as follows:

B i u m i -g i = (B pde i u m i -g pde i ) + (B ic i u m i -g ic i ), (64)
where we have written B i = B pde i + B ic i and g i = g pde i + g ic i to distinguish the contribution of the differential operator from that of the initial condition. Our results (not displayed for brevity) show that after a few greedy iterations, the two contributions have about the same size. Moreover, as the greedy iteration approaches convergence, there is some compensation between the two contributions to the relative residual in Method 1 (but not for Method 2) since they have a size which is about one order of magnitude larger than the relative residual itself. As a further comparison, we considered the quantities

r m i = B 1 u m i -g 1 2 / g 1 2 , ( 65
)
which represent the relative residual for the linear system originating from Method 1 when the iterates produced by Method i ∈ {1, 2, 3} are inserted into the residual. As expected from the MinRes formulation, r m 1 ≤ min(r m 2 , r m 3 ) for all m ≥ 0, and as the greedy iterations approach convergence, r m 1 reaches the value 4 × 10 -5 , whereas r m 2 and r m 3 reach a value of 4 × 10 -4 and 10 -4 , respectively. Figure 2 presents the first six space and time modes for Method 1. The first six modes obtained with Method 2 are essentially the same, whereas some differences, especially in the space modes, can be observed with Method 3. This indicates that the preconditioner plays a relevant role in the exploration of the discrete trial space. Figure 3 reports the normalized errors (u m hk,i -u) as a function of the iteration counter m of the greedy algorithm where, as above, the additional subscript i ∈ {1, 2, 3} indicates which method has been used. The errors are measured in the L 2 (I; H 1 (Ω))-and H 1 (I; H -1 (Ω))-norms. The three methods produce in both norms relatively close errors, and the error for Method 3 is always a bit larger. At convergence of the greedy algorithm, both errors are very small, namely 2 × 10 -4 in the L 2 (I; H 1 (Ω))-norm and 5 × 10 -4 in the H 1 (I; H -1 (Ω))norm.

Figure 4 presents a convergence study for Methods 1, 2 and 3 as a function of the discretization parameters N h (in space) and N k (in time). In both cases, we report the L 2 (I; H 1 (Ω))-and H 1 (I; H -1 (Ω))-errors. The left panel considers N h = (2 l ) 2 , l ∈ {2, 3, 4} with N k = 2 13 , whereas the right panel considers N k = 2 l , l ∈ {4, . . . , 11} with N h = (2 6 ) 2 . We observe that the convergence is of second order in time (if the time-step is small enough) and first order in space in the L 2 (I; H 1 (Ω))-norm, whereas it is of first order in time (if the time-step is small enough) and second order in space in the H 1 (I; H -1 (Ω))-norm. These convergence orders are consistent with the expected decay rates of the bestapproximation errors in both norms when approximating smooth functions by elements of the discrete trial space X hk . Incidentally, we observe that the errors produced by Method 2 in both norms are slightly worse for the coarser time discretizations; this observation is consistent with the CFL-dependent inf-sup stability estimate for Method 2.

Test case 2: time-dependent diffusion

We consider a time-dependent, selfadjoint differential operator A(t) = -µ(t)∆ with diffusion coefficient µ(t) = sin(100πt) + 2. The initial condition is u 0 = 0 and the source term is f = 1. The explicit expression of the exact solution is not available. The discretization parameters are N h = (2 6 ) 2 and N k = 2 13 (as in the previous test case), and the stopping tolerances are greedy = 10 -5 and alt = 5 × 10 -2 (as in the previous test case). The left panel of Figure 5 presents the decrease of the relative residual as a function of the number of greedy iterations for Methods 1, 2 and 3. We notice that the greedy algorithm takes between 16 and 21 iterations for the three methods to converge. The right panel of Figure 5 presents the cumulated number of alternating minimization iterations in the greedy algorithm for Methods 1, 2 and 3. We observe that this number is about the same for Methods 1 and 2 (as for test case 1), whereas it is about 1.5 times larger for Method 3, confirming once again the benefit of using a preconditioner. It is interesting to notice that with Methods 1 and 2, we have solved at convergence of the greedy algorithm about 80 linear systems in space, which is about 0.5% of the amount that would have been solved by using an implicit time-stepping method (recall that N k = 2 13 ). Furthermore, similar observations as for test cases 1 and 2 can be made concerning the two contributions to the relative residual and the behavior of the residuals r m i defined by (65). In particular, we have again r m 1 ≤ min(r m 2 , r m 3 ) for all m ≥ 0 (as expected from the MinRes formulation); as the greedy algorithm approaches convergence, r m 1 reaches a value of 9 × 10 -5 , whereas r m 2 and r m 3 reach a value of 10 -4 and 2 × 10 -4 , respectively. Figure 6 presents the first six space and time modes for Method 1. The first six modes obtained with Method 2 are essentially the same, whereas some differences, especially in the space modes, can be observed with Method 3. This observation again confirms that the preconditioner plays a relevant role in the exploration of the discrete trial space.

Figure 7 reports the normalized differences (u m hk,1 -u m hk,2 ) and (u m hk,1 -u m hk,3 ) as a function of the iteration counter m of the greedy algorithm, where, as above, the additional subscript i ∈ {1, 2, 3} indicates which method has been used. These differences are measured in the L 2 (I; H 1 (Ω))-and H 1 (I; H -1 (Ω))norms. We observe that the three methods produce approximate solutions that are relatively close in both norms. At the convergence of the greedy algorithm, the difference in the L 2 (I; H (Ω))-norm is of the order of 5 × 10 -5 , and it is about one order of magnitude higher in H 1 (I; H -1 (Ω))-norm. Figure 8 presents convergence study for Methods 1, 2 and 3 as a function of the discretization parameters N (in space) and N k (in time). In both cases, we report the L 2 (I; H 1 (Ω))-and H 1 (I; H -1 (Ω))-errors. The left panel considers 13 , the right panel N k = 2 l , l ∈ {4, . . . , 10} with N h = (2 6 ) 2 . Since the exact solution is not available, we consider for each method the approximate solution produced on the finest space-time discretization available. These method-dependent reference solutions are very close according to Figure 7, and their differences are, in both norms, two orders of magnitude lower than the convergence errors reported in Figure 8. In this figure, we observe that for the three methods, the convergence rates are consistent with the best-approximation properties of the discrete trial space X hk in both norms: the convergence is of second order in time and first order in space in the L 2 (I; H 1 (Ω))-norm, whereas it is of first order in time and second order in space in the H 1 (I; H -1 (Ω))-norm.

N h = (2 l ) 2 , l {2, 3, 4} with N k = 2

Test case 3: advection-diffusion

We consider in this section a time-independent, but non-selfadjoint, differential operator A = -∇•µ∇+c(x, y)•∇ with diffusion coefficient µ = 0.1 and advection velocity field c(x, y) = 2π( 12 -y, x-1 2 ) T . The source term is f = 0 and the initial condition is u 0 (x, y) = exp -

(x-2 3 ) 2 +(y-1 2 ) 2 0.07 2
. The explicit expression of the exact solution is not available. The discretization parameters are N h = (2 5 ) 2 and N k = 2 10 , and the stopping tolerances are greedy = 10 -5 and alt = 5×10 -2 (as in the previous test cases). The discretization parameters for this test case are a bit coarser than for the two other test cases because this test case turns out to be more computationally intensive. This is due to the fact that the differential operator in space is non-selfadjoint so that it is necessary to assemble the matrix A T h D -1 h appearing in the definition (38) of the global system matrix B (recall that P = 1 here and that A T h = D h when the differential operator corresponds to a pure diffusion operator).

The left panel of Figure 9 presents the decrease of the relative residual as a function of the number of greedy iterations for Methods 1, 2 and 3. We notice that for Methods 1 and 2, the greedy algorithm takes around 90 iterations to converge (94 and 97, respectively), whereas it takes only 49 iterations for Method 3. Thus, for this test case, Method 3 takes less iterations. The right panel of Figure 9 presents the cumulated number of alternating minimization iterations in the greedy algorithm for Methods 1, 2 and 3. We observe that this number is about the same for the three methods. When reaching convergence for the greedy algorithm, we have solved about 750 linear systems in space, which is 73% of the amount that would have been solved by using an implicit timestepping method (recall that N k = 2 10 ). This percentage is larger than the ones reported for the previous two test cases, but is still competitive. Furthermore, similar observations as for test cases 1 and 2 can be made concerning the two contributions to the relative residual and the behavior of the residuals r m i defined by (65). In particular, we have again r m 1 ≤ min(r m 2 , r m 3 ) for all m ≥ 0 (as expected from the MinRes formulation); as the greedy algorithm approaches convergence, r m 1 reaches a value of 2 × 10 -5 , whereas r m 2 and r m 3 reach a value of 7 × 10 -5 and 3 × 10 -3 , respectively.

Figure 10 presents the first six space and time modes for Method 1. The first six modes obtained with Method 2 are essentially the same, whereas some (left column) and various time-steps N -1 k (right column); in all cases, the curves for Methods 1 and 3 overlap. differences, especially in the space modes, can be observed with Method 3. This observation again confirms that the preconditioner plays a relevant role in the exploration of the discrete trial space.

Figure 11 reports the normalized differences (u m hk,1 -u m hk,2 ) and (u m hk,1 -u m hk,3 ) as a function of the iteration counter m of the greedy algorithm where, as above, the additional subscript i ∈ {1, 2, 3} indicates which method has been used. These differences are measured in the L 2 (I; H 1 (Ω))-and H 1 (I; H -1 (Ω))-norms. We observe that the difference between the solutions produced by Methods 1 and 3 is significant in both norms (two orders of magnitude larger than the difference between Methods 1 and 2); therefore, we can conclude that Method 3 converges more rapidly than Methods 1 and 2, but with a poorer accuracy.

Figure 12 presents a convergence study for Methods 1, 2 and 3 as a function of the discretization parameters N h (in space) and N k (in time). In both cases, we report the L 2 (I; H 1 (Ω))-and H 1 (I; H -1 (Ω))-errors. The left panel considers N h = (2 l ) 2 , l ∈ {2, 3, 4} with N k = 2 10 , whereas the right panel considers N k = 2 l , l ∈ {4, . . . , 8} with N h = (2 5 ) 2 . Since the exact solution is not available, we consider for each method the approximate solution produced on the finest space-time discretization available. For Methods 1 and 2, these two reference solutions are very close according to Figure 11, and their difference is, in both norms, two orders of magnitude lower than the convergence errors reported in Figure 12. Moreover, for both methods, the reported convergence rates are, as above, consistent with the best-approximation properties of the discrete trial space X hk in both norms. Finally, for Method 3, the convergence rates are similar except for the behavior with respect to time refinement in the L 2 (I; H 1 (Ω))-norm which is somewhat sub-optimal.

Conclusion and outlook

In this work, we have devised a space-time tensor for the low-rank approximation of linear parabolic evolution equations. The proposed method is uniformly stable with respect to the space-time discretization parameters and leads to solving sequentially global problems in space and in time. Our numerical results on various test cases indicate the importance of the preconditioner (which enters the method by means of the norm equipping the discrete test space), since the preconditioner has an influence on the convergence of the greedy algorithm. However, although the preconditioner improves the convergence of the greedy algorithm, it increases the cost of each iteration. Table 1 summarizes the relative CPU times for Methods 2 and 3 normalized with respect to Method 1 (we have considered 21 random initializations in each case and used for each method the median CPU time). We can see from this table that Methods 1 and 2 essentially deliver the same CPU times, whereas Method 3 turns out to be more effective especially for test case 3 despite the increased number of greedy iterations. Finally, various perspectives of this work can be envisaged. We mention in particular the question of adapting the discretization spaces and that of devising different approaches to obtain a separated representation of the exact solution with sufficient accuracy at low-rank when the differential operator has

Remark 4 . 1 (

 41 Matrix form). The matrix form of problem (43) is as follows:

Figure 1 :

 1 Figure 1: Test case 1. Left: relative residual at each iteration of the greedy algorithm. Right: cumulated number of alternating minimization iterations in the greedy algorithm.

Figure 2 :

 2 Figure 2: Test case 1: first six modes in space (top row) and in time (bottom row) for Method 1.

Figure 3 :

 3 Figure 3: Test case 1: comparison of the errors produced by Methods 1, 2, 3 in two norms: L 2 (I; H 1 (Ω)) (left) and H 1 (I; H -1 (Ω)) (right); in both cases, the curves for Methods 1 and 2 almost overlap.

Figure 4 : 1 k

 41 Figure 4: Test case 1: convergence study for Methods 1, 2 and 3 for errors measured in the L 2 (I; H 1 (Ω))-norm (top row) and in the H 1 (I; H -1 (Ω))-norm (bottom row) for various mesh-sizes N -1/2 h (left column) and various time-steps N -1 k (right column); in all cases, the curves for Methods 1 and 3 overlap.

Figure 5 :

 5 Figure 5: Test case 2. Left: relative residual at each iteration of the greedy algorithm. Right: cumulated number of alternating minimization iterations in the greedy algorithm.

Figure 6 :

 6 Figure 6: Test case 2: first six modes in space (top row) in time (bottom row) for Method 1.

Figure 7 :

 7 Figure 7: Test case 2: comparison of the solutions produced by Methods 1, 2, 3 in two norms: L 2 (I; H 1 (Ω)) (left) and H 1 (I; H -1 (Ω)) (right).

22

 22

Figure 8 :

 8 Figure 8: Test case 2: convergence study for Methods 1, 2 and 3 for errors measured in the L 2 (I; H 1 (Ω))-norm (top row) and in the H 1 (I; H -1 (Ω))-norm (bottom row) for various mesh-sizes N -1/2 h

Figure 9 :

 9 Figure 9: Test case 3. Left: relative residual at each iteration of the greedy algorithm. Right: cumulated number of alternating minimization iterations in the greedy algorithm.

Figure 10 :

 10 Figure 10: Test case 3: first six modes in space (top row) and in time (bottom row) for Method 1.

Figure 11 :

 11 Figure 11: Test case 3: comparison of the solutions produced by Methods 1, 2, 3 in two norms: L 2 (I; H 1 (Ω)) (left) and H 1 (I; H -1 (Ω)) (right).

Figure 12 : 1 k

 121 Figure 12: Test case 3: convergence study for Methods 1, 2 and 3 for errors measured in the L 2 (I; H 1 (Ω))-norm (top row) and in the H 1 (I; H -1 (Ω))-norm (bottom row) for various mesh-sizes N -1/2 h (left column) and various time-steps N -1 k (right column); in the right column, the curve for Method 2 overlaps with that of Method 1 for the smaller time-steps and with that of Method 3 for the larger time-steps.