Stability properties of a projector-splitting scheme for dynamical low rank approximation of random parabolic equations

We consider the Dynamical Low Rank (DLR) approximation of random parabolic equations and propose a class of fully discrete numerical schemes. Similarly to the continuous DLR approximation, our schemes are shown to satisfy a discrete variational formulation. By exploiting this property, we establish stability of our schemes: we show that our explicit and semi-implicit versions are conditionally stable under a parabolic type CFL condition which does not depend on the smallest singular value of the DLR solution; whereas our implicit scheme is unconditionally stable. Moreover, we show that, in certain cases, the semi-implicit scheme can be unconditionally stable if the randomness in the system is sufficiently small. Furthermore, we show that these schemes can be interpreted as projector-splitting integrators and are strongly related to the scheme proposed by Lubich et al. [BIT Num. Math., 54:171-188, 2014; SIAM J. on Num. Anal., 53:917-941, 2015], to which our stability analysis applies as well. The analysis is supported by numerical results showing the sharpness of the obtained stability conditions.


Introduction
Many physical and engineering applications are modeled by time-dependent partial differential equations (PDEs) with input data often subject to uncertainty due to measurement errors or insufficient knowledge. These uncertainties can be often described by means of probability theory by introducing a set of random variables into the system. In the present work, we consider a random evolutionary PDEu with random initial condition, random forcing term and a random linear elliptic operator L. Many of the numerical methods used to approximate such problems, require evaluating the, possibly expensive, model in many random parameters.
In this regard, the use of reduced order models (e.g. Proper orthogonal decomposition [5,6] or generalized Polynomial chaos expansion [37,39,26,10,33]) is of a high interest. When the dependence of the solution on the random parameters significantly changes in time, the use of time-varying bases is very appealing. In the present work, we consider the dynamical low rank (DLR) approximation (see [34,32,7,28,4,22,23,16]) which allows both the deterministic and stochastic basis functions to evolve in time while exploiting the structure of the differential equation. An extension to tensor differential equations was proposed in [25,31]. The DLR approximation of the solution is of the form where R is the rank of the approximation and is kept fixed in time,ū(t) = E[u(t)] is the mean value of the DLR solution, {U j (t)} R j=1 is a time dependent set of deterministic basis functions, {Y j (t)} R j=1 is a time dependent set of zero mean stochastic basis functions. By suitably projecting the residual of the differential equation one can derive evolution equations for the mean valueū and the deterministic and stochastic modes {U j } R j=1 , {Y j } R j=1 (see [34,23]), which in practice need to be solved numerically. An efficient and stable discretization scheme is therefore of a high interest.
In [34,23], Runge-Kutta methods of different orders were applied directly to the system of evolution equations for the deterministic and stochastic basis functions. In the presence of small singular values in the solution, the system of evolution equations becomes stiff as an inversion of a singular or nearly-singular matrix is required to solve it. Applying standard explicit or implicit Runge-Kutta methods leads to instabilities (see [20]). In this respect, the projector-splitting integrators (proposed in [29,30] and applied in e.g. [12,11]) are very appealing. In [20], the authors showed that when applying the projector-splitting method for matrix differential equations one can bound the error independently of the size of the singular values, under the assumption that f − L maps onto the tangent bundle of the manifold of all R-rank functions up to a small error of magnitude ε. A limitation of their theoretical result, as the authors point out, is that it requires a Lipschitz condition on f − L and is applicable to discretized PDEs only under a severe condition tL 1 where t is the step size and L is the Lipschitz constant, even for implicit schemes. Such condition is, however, not observed in numerical experiments. In [21], the authors proposed projected Runge-Kutta methods, where following a Runge-Kutta integration, the solution first leaves the manifold of R-rank functions by increasing its rank, and then is retracted back to the manifold. Analogous error bounds as in [20] are obtained, also for higher order schemes, under the same ε-approximability condition on f − L and under a restrictive parabolic condition on the time step.
In this work we propose a class of numerical schemes to approximate the evolution equations for the mean, the deterministic basis and the stochastic basis, which can be of explicit, semi-implicit or implicit type. Although not evident at first sight, we show that the explicit version of our scheme can be reinterpreted as a projector-splitting scheme, whenever the discrete solution is full-rank, and is thus equivalent to the scheme from [29,30]. However, our derivation allows for an easy construction of implicit or semi-implicit versions.
The main goal of this work is to prove the stability of the proposed numerical schemes for a parabolic problem (2). We first show that the continuous DLR solution satisfies analogous stability properties as the weak solution of the parabolic problem (1). We then analyze the stability of the fully discrete schemes. Quite surprisingly, the stability properties of both the discrete and the continuous DLR solutions do not depend on the size of their singular values, even without any ε-approximability condition on f − L. The implicit scheme is proven to be unconditionally stable. This improves the stability result which could be drawn from the error estimates derived in [20]. The explicit scheme remains stable under a standard parabolic stability condition between time and space discretization parameters for an explicit propagation of parabolic equations. The semi-implicit scheme is generally only conditionally stable under again a parabolic stability condition, and becomes unconditionally stable under some restrictions on the size of the randomness of the operator. As an application of the general theory developed in this paper, we consider the case of a heat equation with a random diffusion coefficient. We dedicate a section to particularize the numerical schemes and the corresponding stability results to this problem. The semi-implicit scheme turns out to be always unconditionally stable if the diffusion coefficient depends affinely on the random variables. The sharpness of the obtained stability conditions on the time step and spatial discretization is supported by the numerical results provided in the last section.
A big part of the paper is dedicated to proving a variational formulation of the discretized DLR problem, analogous to the variational formulation of the continuous DLR problem (see [32,Prop. 3.4]). Such formulation is a key for showing the stability properties and, as we believe, might be useful for some further analysis of the proposed discretization schemes. It as well applies to the projector-splitting integrator from [29,30] provided the solution remains full rank at all time steps. However, in the rank-deficient case, our schemes may result in different solutions. We dedicate a subsection to show that a rank-deficient solution obtained by our scheme still satisfies a suitable discrete variational formulation and consequently has the same stability properties as the full-rank case.
The outline of the paper is the following: in Section 2 we introduce the problem and basic notation; in Section 3 we describe the DLR approximation and recall its geometrical interpretation with variational formulation. In Section 4 we describe the discretization of the DLR method and propose three types of time integration schemes. We then derive a variational formulation for the discrete DLR solution and show its reinterpretation as a projector-splitting scheme. Section 5 is dedicated to proving the stability properties of both continuous and discrete DLR solution. In Section 6, we analyze the case of a heat equation with random diffusion coefficient and random initial condition. Finally in Section 7 we present several numerical tests that support the derived theory. Section 8 draws some conclusions.

Problem statement
We start by introducing some notation. Let (Γ, F, ρ) be a probability space. Consider the Hilbert space L 2 ρ = L 2 ρ (Γ) of real valued random variables on Γ with bounded second moments, with associated scalar product v, w L 2 ρ = Γ vw dρ and norm v L 2 ρ = v, v L 2 ρ . Consider as well two separable Hilbert spaces H and V with scalar products ·, · H , ·, · V , respectively. Suppose that H and V form a Gelfand triple (V, H, V ), i.e. V is a dense subspace of H and the embedding V → H is continuous with a continuity constant C P > 0. Let ) is a Gelfand triple as well (see e.g. [27,Th. 8.17]), and we have We define the mean value of a random variable v as E[v] = Γ v dρ, where the integral here denotes the Bochner integral in a suitable sense, depending on the co-domain of the random variable considered. In what follows, we will use the notationv = E[v] and v * := v −v. Moreover, we let (·, ·) V V,L 2 ρ denote the dual pairing between L 2 ρ (Γ; V ) and L 2 ρ (Γ; V ): With this notation at hand, we now consider a random operator L with values in the space of linear bounded operators from V to V that is uniformly bounded and coercive, i.e. a Borel measurable function Associated to the random operator L, we introduce the operator L, defined as Notice that for any strongly measurable u : Γ → V the map ω ∈ Γ → L(ω)u(ω) ∈ V is strongly measurable, V being separable, see Proposition A in the appendix. From the uniform boundedness of L it follows immediately that, if u is square integrable, then L(u) is square integrable as well and L(u) which is coercive and bounded with coercivity and continuity constant C L and Then, given a final time T > 0, a random forcing term f ∈ L 2 (0, T ; L 2 ρ (Γ; H)) and a random initial condition u 0 ∈ L 2 ρ (Γ; V ), we consider now the following parabolic problem: Find a solution u true ∈ L 2 (0, T ; The general theory of parabolic equations (see e.g. [38]) can be applied to problem (6), at least in the case of L 2 ρ (Γ; V ), L 2 ρ (Γ; H), L 2 ρ (Γ; V ) being separable, e.g. when Γ is a Polish space and F is the corresponding Borel σ-algebra. We conclude then that problem (6) has a unique solution u true which depends continuously on f and u 0 . We note that the theory of parabolic equations would allow for less regular data f ∈ L 2 (0, T ; L 2 ρ (Γ; V )) and u 0 ∈ L 2 ρ (Γ; H). However, in this work we restrict our attention to the case f ∈ L 2 (0, T ; L 2 ρ (Γ; H)), u 0 ∈ L 2 ρ (Γ; V ).

Dynamical low rank approximation and its variational formulation
Dynamical low rank (DLR) approximation, or dynamically orthogonal (DO) approximation (see e.g. [23,34,24]) seeks an approximation of the solution u true of problem (6) in the form whereū(t) ∈ V , {U j (t)} R j=1 ⊂ V is a time dependent set of linearly independent deterministic basis functions, {Y j (t)} R j=1 ⊂ L 2 ρ is a time dependent set of linearly independent stochastic basis functions. In what follows, we focus on the so called Dual DO formulation (see e.g. [32]), in which the stochastic basis are only required to be linearly independent at all times. We call R the rank of a function u of the form (7). To ensure the uniqueness of the expansion (7) for a given initialization u(0) =ū(0) + R j=1 U j (0)Y j (0), we consider the following conditions: and the gauge condition (also called DO condition) (see [19]). Plugging the DLR expansion (7) into the equation (6) and following analogous steps as proposed in [34] leads to the DLR system of equations presented next.
Definition 3.1 (DLR solution). We define the DLR solution of problem (6) as are solutions of the following system of equations: with the initial conditionsū(0), is a good approximation of u 0 . In (12), the matrix M ∈ R R×R is defined as M ij := U i , U j H , 1 ≤ i, j ≤ R and P ⊥ Y denotes the orthogonal projection operator in the space L 2 ρ (Γ) on the orthogonal complement of the R-dimensional subspace Y = span{Y 1 , . . . , Y R }, i.e.
For the initial condition one can use for instance a truncated Karhunen- are the first R (rescaled) eigenfunctions of the covariance operator C u0 : H → H defined as In what follows we will use the notation U = (U 1 , . . . , U R ) and Y = (Y 1 , . . . , Y R ). Then, the approximation (7) reads u =ū + U Y .
The rest of the section gives a geometrical interpretation of the DLR method and derives a variational formulation, following to a large extent derivations from [32]. Such geometrical interpretation and consequent variational formulation will be key to derive the stability results of the numerical schemes, discussed in Section 5.2. We first introduce the notion of a manifold of R-rank functions, characterize its tangent space in a point as well as the orthogonal projection onto the tangent space.
The vector space consisting of all square integrable random variables with zero mean value will be denoted by L 2 ρ,0 = L 2 ρ,0 (Γ) ⊂ L 2 ρ (Γ). Definition 3.2 (Manifold of R-rank functions). By M R ⊂ L 2 ρ,0 (Γ; V ) we denote the manifold consisting of all rank R random functions with zero mean It is well known that M R admits an infinite dimensional Riemannian manifold structure ( [15]).
where U = span{U 1 , . . . , U R } and P U [·] is the H-orthogonal projection onto the subspace U.
For more details, see e.g. [32]. Note that Π U Y [·] can be equivalently written . In the following we will extend the domain of the projection operator Π U Y . Further, we will state two lemmas used to establish Theorem 3.7, which presents the variational formulation of the DLR approximation.
The operator Π U Y can be extended to an operator from L 2 . The extended operator satisfies the following.
Proof. First, we show that Indeed, where in the forth step we applied Theorem 8.13 from [27]. Now we proceed with proving (17) We are now in the position to state the first variational formulation of the DLR equations. Lemma 3.6. Let U, Y be the solution of the system (11)- (12). Then the zeromean part of the DLR solution u * = U Y satisfies Proof. First, we multiply equation (11) by Y j and take its weak formulation in L 2 ρ . Summing over j results in Analogously, we multiply (12) by U j and take its weak formulation in V Summing over j, this leads to Summing the derived equations we obtain In particular, this holds for any z being a Bochner integrable simple function, the collection of which is dense in L 2 ρ (Γ; V ) (see [27,Th. 8.15]).
We can finally state the variational formulation corresponding to the DLR equations (10)- (12).
Theorem 3.7 (DLR variational formulation). Letū, U, Y be the solution of the system (10)- (12). Then the DLR solution u =ū + U Y satisfies Proof. Based on Lemma 3.6 and Lemma 3.5 we can write which can be equivalently written as exploiting the fact that u * + L * (u) − f * , w V V,L 2 ρ = 0, ∀w ∈ V . Likewise, equation (10) can be equivalently written as Summing (20) and (21) leads to the sought equation (19).
Recently, the existence and uniqueness of the dynamical low rank approximation for a class of random semi-linear evolutionary equations was established in [19] and for linear parabolic equations in two space dimensions with a symmetric operator L in [3].

Discretization of DLR equations
In this section we describe the discretization of the DLR equations that we consider in this work. In particular, we focus on the time discretization of (10)- (12) and propose a staggered time marching scheme that decouples the update of the spatial and stochastic modes. Afterwards, we will show that the proposed scheme can be formulated as a projector-splitting scheme for the Dual DO formulation and comment on its connection to the projector-splitting scheme from [29]. As a last result we state and prove a variational formulation of the discretized problem.
Stochastic discretization We consider a discrete measure given by {ω k , λ k }N k=1 , i.e. a set of sample points {ω k }N k=1 ⊂ Γ with R <N < ∞ and a set of positive weights {λ k }N k=1 , λ k > 0, N k=1 λ k = 1, which approximates the probability measure ρρ The discrete probability space (Γ = {ω k }N k=1 , 2Γ,ρ) will replace the original one (Γ, F, ρ) in the discretization of the DLR equations. Notice, in particular, that a random variable Z :Γ → R measurable on (Γ, 2Γ,ρ) can be represented as a vector z ∈ RN with z k = Z(ω k ), k = 1, . . . ,N . The sample points {ω k }N k=1 can be taken as iid samples from ρ (e.g. Monte Carlo samples) or chosen deterministically (e.g. deterministic quadrature points with positive quadrature weights). The mean value of a random variable Z with respect to the measurê ρ is computed as We introduce also the semi-discrete scalar products ·, · ,L 2 ρ with = V, H and their corresponding induced norms · ,L 2 ρ . Note that the semi-discrete bilinear form ·, · L,ρ defined as is coercive and bounded, with the same coercivity and continuity constants C L , C B , defined in (4), (5), respectively.
Space discretization We consider a general finite-dimensional subspace V h ⊂ V whose dimension is larger than R and is determined by the discretization parameter h. Eventually, we will perform a Galerkin projection of the DLR equations onto the subspace V h . We further assume that an inverse inequality of the type holds for some p ∈ N and C I > 0.
Time discretization For the time discretization we divide the time interval into N equally spaced subintervals 0 = t 0 < t 1 < · · · < t N = T and denote the time step by t := t n+1 − t n . Note that the DLR solution u =ū + U Y appears in the right hand side of the system of equations (10)- (12) both in the operator L and in the projector operator onto the tangent space to the manifold. We will treat these two terms differently. Concerning the projection operator, we adopt a staggered strategy, where, given the approximate solution u n =ū n + U n Y n , we first update the meanū n+1 , then we update the deterministic basis U n+1 projecting on the subspace span{Y n }; finally, we update the stochastic basis Y n+1 projecting on the orthogonal complement of span{Y n } and on the updated subspace span{U n+1 }. This staggered strategy resembles the projection splitting operator proposed in [29]. We will show later in Section 4.4 that it does actually coincide with the algorithm in [29]. Concerning the operator L, we will discuss hereafter different discretization choices leading to explicit, semi-implicit or fully implicit algorithms.

Fully discrete problem
We give in the next algorithm the general form of the discretization schemes that we consider in this work.
1. Compute the mean valueū n+1 such that is the analogue of the projector defined in (13) but in the discrete space L 2 ρ .

Reorthonormalize the stochastic basis: find
5. Form the approximated solution at time step t n+1 as The expressions L(u n h,ρ , u n+1 h,ρ ) and f n,n+1 stand for an unspecified time integration of the operator L(u(t)) and right hand side f (t), t ∈ [t n , t n+1 ] and v * denotes the 0-mean part of a random variable v ∈ L 2 ρ with respect to the discrete The newly computed solution u n+1 h,ρ belongs to the tensor product space (25) can be rewritten as a deterministic linear system of R×N equations with R×N unknowns. This system can be decoupled into a linear system of size R × R for each collocation point. If the deterministic modesŨ n+1 are linearly independent, the system matrix is invertible. Otherwise we interpret (25) in a minimal-norm least squares sense, choosing a solutioñ Y n+1 , if it exists, that minimizes the norm Ỹ n+1 − Y n L 2 ρ . This is discussed in more details in Section 4.3.
The following lemma shows that the scheme (23)-(25) satisfies some important properties that will be essential in the stability analysis presented in Section 5.

Lemma 4.2 (Discretization properties).
Assuming that a solution (Ỹ n+1 ,Ũ n+1 ,ū n+1 ) exists, the following properties hold for the discretization (23)-(25): 1. Discrete DO condition: 1. In the following proof we assume that the matrixM n+1 = Ũ n+1 ,Ũ n+1 H is full rank. For the rank-deficient case, we refer the reader to the proof of Lemma 4.12. Let us multiply equation (25) by Y n from the left and take the L 2 ρ -scalar product. Since the second term involves P ⊥ ρ,Y n , the scalar product of Y n with the second term vanishes which, under the assumption thatM n+1 is full rank, gives us the discrete DO condition 2. This is a consequence of the fact that we have Eρ[Y n ] = 0 and 3. This is immediate from the discrete DO property and Y n , Y n To complete the discretization scheme (23)-(25) we need to specify the terms L(u n h,ρ , u n+1 h,ρ ) and f n,n+1 . The DLR system stated in (10)-(12) is coupled. Therefore, an important feature we would like to attain is to decouple the equations for the mean value, the deterministic and the stochastic modes as much as possible. We describe hereafter 3 strategies for the discretization of the operator evaluation term L(u n h,ρ , u n+1 h,ρ ), and the right hand side f n,n+1 .
Explicit Euler scheme The explicit Euler scheme performs the discretization It decouples the system (23)-(25) since, for the computation of the new modes, we require only the knowledge of the already-computed modes. The equations for the stochastic modes {Ỹ n+1 j } R j=1 are coupled together through the matrix M n+1 = Ũ n+1 ,Ũ n+1 H ∈ R R×R but are otherwise decoupled between collocation points (i.e.N linear systems of size R have to be solved).
Implicit Euler scheme The implicit Euler scheme performs the discretization This method couples the system (23)-(25) in a non-trivial way, which is why we do not focus on this method in our numerical results. We mention it in the stability estimates section (Section 5.2) for its interesting stability properties.
Semi-implicit scheme Assume that our operator L can be decomposed into two parts where L det : V → V is a linear deterministic operator such that it induces a bounded and coercive bilinear form ·, and that its action on a function Then, L det is also a linear operator L det :

and induces a bounded coercive bilinear form on
We propose a semi-implicit time integration of the operator evaluation term whereas for f n,n+1 we can either take f n,n+1 = f (t n+1 ) or f n,n+1 = f (t n ) or any convex combination of both. The resulting scheme is detailed in the next lemma.
Proof. The equation for the mean (23) using the semi-implicit scheme (30) can be written as Noticing that gives us equation (31). Concerning the equation for the deterministic modes we derive The term T 3 vanishes since Eρ[Y n ] = 0 and the term T 4 can be further expressed as where we used the discrete DO condition (28). Finally, the stochastic equation (25) can be written as The term T 5 vanishes since L * det (ū n+1 ) = 0. As for T 6 , we derive which leads us to the sought equation (33).
We see from (31)-(33) that, similarly to the explicit Euler scheme, the equations for the mean, deterministic modes and stochastic modes are decoupled. If the spatial discretization of the PDEs (31) and (32) is performed by the Galerkin approximation, the final linear system involves the inversion of the matrix where {ϕ i } is the basis of V h in which the solution is represented. Both the mass matrix ϕ j , ϕ i H and the stiffness matrix ϕ j , ϕ i L det are positive definite and do not evolve with time, so that an LU factorization can be computed once and for all at the beginning of the simulation. Concerning the stochastic equation (33), we need to solve a linear system with the matrixM n+1 + t Ũ n+1 ,Ũ n+1 L det for each collocation point ω k , unlike the explicit Euler method, where the system involves only the matrixM n+1 . The matrixM n+1 + t Ũ n+1 ,Ũ n+1 L det is symmetric and positive definite with the smallest singular value bigger than that ofM n+1 . Notice, however, that ifM n+1 is rank deficient, also the matrix M n+1 + t Ũ n+1 ,Ũ n+1 L det will be so. Note that there exists a unique discrete DLR solution for the explicit and semi-implicit version of Algorithm 4.1 also in the rank-deficient case (see Lemma 4.13 below). The existence of solutions for the implicit version remains still an open question.

Discrete variational formulation for the full-rank case
This subsection will closely follow the geometrical interpretation introduced in Section 3. We will introduce analogous geometrical concepts for the discrete setting, i.e. manifold of R-rank functions, tangent space and orthogonal projection, and will show in Theorem 4.9 that the scheme from Algorithm 4.1 can be written in a (discrete) variational formulation, assuming that the matrixM n+1 stays full-rank.
we denote the manifold of all rank R functions with zero mean that belong to the (possibly finite dimensional) space The projection Π h,ρ U Y is defined in the discrete space V h ⊗ L 2 ρ analogously to its continuous version (16). It holds A discrete analogue of Lemma 3.5 holds, i.e.
The solution of the proposed numerical scheme (23)-(26) satisfies a discrete variational formulation analogous to the variational formulation (19). To show this, we first present a technical lemma which will be important in deriving the variational formulation as well as in the stability analysis presented in Section 5.
Lemma 4.6. Let u n h,ρ , u n+1 h,ρ be the discrete DLR solution at t n , t n+1 , respectively, from the scheme in Algorithm 4.1. Then the zero-mean parts u n, * h,ρ , u n+1, * Proof.
1. The solution u n, * h,ρ can be written as Since 0 , Y n L 2 ρ = 0, using the definition (35) we have 2. The newly computed solution u n+1, * h,ρ can be expressed as Based on (28) Remark 4.7. Note that for any function of the form R is a vector space, it includes any linear combination of u n, * h,ρ and u n+1, * h,ρ . The following lemma is an analogue of Lemma 3.6 and will become useful when we derive the discrete variational formulation.
Lemma 4.8. Let u n h,ρ , u n+1 h,ρ be the discrete DLR solutions at times t n , t n+1 as defined in Algorithm 4.1. Then the zero-mean parts u n+1 * h,ρ , u n * h,ρ satisfy Proof. Multiplying (24) by Y n j and summing over j, we obtain Noticing that and taking the weak formulation of (38) Similarly, multiplying (25) byŨ n+1 , and further writing (25) in a weak form in L 2 ρ , we obtain Since taking the weak formulation of (40) Finally, summing equations (39) and (41) results in (37).
We now proceed with the discrete variational formulation.
Theorem 4.9 (Discrete variational formulation). Let u n h,ρ and u n+1 h,ρ be the discrete DLR solution at times t n , t n+1 , respectively, n = 0, . . . , N − 1, as defined in Algorithm 4.1. Then it holds Proof. Thanks to Lemma 4.6 we have (u n+1 h,ρ − u n h,ρ ) * ∈ TŨ n+1 Y n M h,ρ R , and we can derive and formula (36) gives us Summing (43), (44) and applying Lemma 4.8 results in which is equivalent to the final result (42). In (45) we have employed The preceding theorem applies to a discretization of any kind of the operator L ∈ L 2 ρ (Γ; V ), not necessarily elliptic or linear, as assumed in Section 2, as long as Lemma 3.5 holds.

Discrete variational formulation for the rank-deficient case
The discrete variational formulation established in the previous section is valid only in the case of the deterministic basisŨ n+1 being linearly independent, since the proof of Theorem 4.9 implicitly involves the inverse ofM n+1 = Ũ n+1 ,Ũ n+1 H . In this subsection, we show that a discrete variational formulation can be generalized for the rank-deficient case.
When applying the discretization scheme proposed in step 3 of Algorithm 4.1 with a rank-deficient matrixM n+1 , we recall that the solutionỸ n+1 is defined as the solution of (25) minimizing Ỹ n+1 − Y n L 2 ρ . Note that minimizing Ỹ n+1 − Y n L 2 ρ is equivalent to minimizing the norm Ỹ n+1 (ω k ) − Y n (ω k ) R R for every sample point ω k , k = 1, . . . ,N , where · 2 R R = ·, · R R denotes the Euclidean scalar product in R R .
In what follows we will exploit the fact that the vector space L 2 ρ is isomorphic to RN . In particular, it holds that (Ỹ n+1 − Y n ) ∈ R R×N , where each column of (Ỹ n+1 − Y n ) is given by (Ỹ n+1 − Y n )(ω k ), k = 1, . . . ,N . With a little abuse of notation, we useŨ n+1 : R R → V h to denote a linear operator which takes real coefficients and returns the corresponding linear combination of the basis functionsŨ n+1 . ByŨ n+1 : V h → R R we denote its dual.
Proof. Seeking a contradiction, let us suppose that where P ker(M n+1 ) [v] ∈ R R×N for v ∈ R R×N denotes the column-wise application of ·, · R R -orthogonal projection onto the kernel ofM n+1 . Then, such constructed Z satisfies and solves (25): where in the last step we used that ker(M n+1 ) = ker(Ũ n+1 ). This leads to a contradiction thatỸ n+1 was the solution minimizing Ỹ n+1 − Y n L 2 ρ . When showing the equivalence between the DLR variational formulation (19) and the DLR system of equations (10)- (12) in the continuous setting, the DO condition (9) plays an important role. In an analogous way, the discrete DO condition (property 1 from Lemma 4.2 for the full-rank case) plays an important role when showing the equivalence between the discrete DLR system of equations and the discrete DLR variational formulation.
Proof. LetỸ n+1 be a solution of (25) minimizing Ỹ n+1 − Y n L 2 ρ . Thanks to Lemma 4.11 we know that for any v ∈ ker(M n+1 ) ⊥ , the solutionỸ n+1 of equation (25) satisfies then the statement will follow. But for the column space of ⊂ RN being the orthogonal complement to Y n in the scalar product ·, · L 2 ρ . Now the proof is complete. In the following lemma we address the question of existence of a unique solution when applying the explicit and semi-implicit scheme. Proof. We will start with the semi-implicit scheme. By virtue of Lemma 4.3, under the discrete DO condition (47), applying the semi-implicit scheme to equation (25) is equivalent to solving equation (33). We will first focus our attention to equation (33) and show that there exists a unique solution minimizing Ỹ n+1 − Y n L 2 ρ . This solution will satisfy the discrete DO and consequently is a unique minimizing solution of (25). Equation (33) can be rewritten as where Since RHS above lies in the range ofŨ n+1 , which is the same as the range of B, a solution of (49) exists. Moreover, since the matrix B is positive definite on the space ker(B) ⊥ , any solution can be expressed as (Ỹ n+1 − Y n + W ) with W ∈ ker(B) N and a uniqueỸ n+1 ∈ R R×N such that (Ỹ n+1 − Y n ) ∈ ker(B) ⊥ N . The solutionỸ n+1 minimizes each column (Ỹ n+1 − Y n )(ω k ) R R , k = 1, . . . ,N and thus it is the unique solution of (49) that minimizes norm Ỹ n+1 − Y n L 2 ρ . We observe that the established solutionỸ n+1 of equation (49) satisfies the discrete DO condition (47). The argument is analogous to the proof of Lemma 4.12, but instead ofM n+1 here we take B. Therefore, the statement for the semi-implicit scheme follows. The explicit case can be shown by following analogous steps with Now we can proceed with showing the discrete variational formulation. It is not generally easy to deal with the notion of a tangent space at a certain point on the manifold in the rank-deficient case. In the following theorem we will, however, show that an analogous discrete variational formulation holds. Given U ∈ (V h ) R and Y ∈ (L 2 ρ,0 ) R , we define the vector space T U Y as It is easy to verify that, analogously to Lemma 4.6, the (possibly rank-deficient) discrete DLR solutions u n h,ρ and u n+1 h,ρ at times t n , t n+1 , as defined in Algorithm 4.1 satisfy Theorem 4.14. Let u n h,ρ and u n+1 h,ρ be the (possibly rank-deficient) discrete DLR solution at times t n , t n+1 , respectively, n = 0, . . . , N − 1, as defined in Algorithm 4.1. Then the following variational formulation holds Proof. First, consider equation (24) with v h =Ũ n+1 j . Summing over j results in (52) Let us proceed with the equation (25): Taking a weak formulation in L 2 ρ,0 results in Concerning equation (24), we proceed as follows: where in the second step we applied Eρ[Ỹ n+1 Y n ] = Id which holds thanks to the discrete DO condition from Lemma 4.12. Summing equation (53) and (54) we obtain The rest of the proof follows the same steps as in the proof of Theorem 4.9, i.e. summing the mean value equation (23) and noting that some terms vanish.
Remark 4.15. Thanks to the observation that ker(M n+1 ) = ker(Ũ n+1 ), we can easily see that any discrete solutionỸ n+1 of equation (25) leads to the same discrete DLR solution u n+1 h,ρ =ū n+1 h,ρ +Ũ n+1Ỹ n+1 . Therefore, the result of the preceding theorem as well as the stability properties shown in Section 5 hold for the discrete DLR solution obtained by any of the solutions of equation (25).

Reinterpretation as a projector-splitting scheme
The proposed Algorithm 4.1 was derived from the DLR system of equations (10)- (12). This subsection is dedicated to showing that this scheme can in fact be formulated as a projector-splitting scheme for the time discretization of the Dual DO approximation of (6). Afterwards, we will continue by showing its connection to the projector-splitting scheme of the first order proposed in [29,30] and further analyzed in [20].
In what follows, we will focus on the evolution of u n, * h,ρ , i.e. the 0-mean part of the discrete DLR solution u n h,ρ .
Lemma 4.16. The discretized system of equations (24)-(25) can be equivalently reformulated as Proof. These equations are essentially equations (39) and (41), which are shown to hold in the proof of Lemma 4.8.
We recall that from Lemma 3.6, the zero-mean part of the continuous DLR approximation u * = U Y satisfies

Comparison to the projection scheme in [29]
There are several equivalent DLR formulations. The DO formulation, proposed and applied in [34,35,36], seeks for an approximation of the form u R = U Y with {U j } R j=1 ⊂ V h orthonormal in ·, · H and {Y j } R j=1 ⊂ L 2 ρ linearly independent. The dual DO formulation, on the contrary, keeps the stochastic basis {Y j } R j=1 orthonormal in ·, · L 2 ρ and {U j } R j=1 linearly independent. The double dynamically orthogonal (DDO) or bi-orthogonal formulation searches for an approximation in the form orthonormal in ·, · H , ·, · L 2 ρ , respectively, and S ∈ R R×R a full rank matrix (see e.g. [7,8,23]). In [9,32] it was shown that these formulations are equivalent. In our work we consider the dual DO formulation with an isolated mean so that the stochastic basis functions are centered.
A first order projector-splitting scheme introduced in [29,30] and further analyzed in [20] is a time integration scheme successfully used for the integration of dynamical low rank approximation in the DDO formulation. This subsection provides a detailed look into the comparison of the Algorithm 4.1 and the discretization scheme from [29,30]. We will see that, if the solution is full rank, these schemes are in fact equivalent.
We will adapt the algorithm from [29] to approximate the DLR solution in the DDO form with an isolated mean, i.e.
Having an R-rank solution u n h,ρ , the basic first-order scheme from [29] requires the knowledge of the solution u n+1 h,ρ , which is used in evaluating the term A = u n+1 h,ρ − u n h,ρ . To deal with differential equations where u n+1 h,ρ is a-priori unknown, we will consider a general scheme where where f n,n+1 and L(u n h,ρ , u n+1 h,ρ ) can be any of the explicit, implicit or semiimplicit discretizations detailed in Section 4.1. Adopting the notation from [29], the splitting scheme from [29,30] for a DDO approximation of (6) results in the following 6-step algorithm. 1. Compute the mean valueû n+1 such that 2. Solve for K 1 such that

SetS
The new solutionû n+1 h,ρ is then defined aŝ Now, let us compare the previous steps to Algorithm 4.1. We can easily observe thatû n+1 =ū n+1 . Since Y n = V 0 , we can see that equation (24) is equivalent to step 1 with U n = U 0 S 0 , i.e. K 1 =Ũ n+1 . Further, we havẽ Equation (25) can be reformulated as Note that the expression in brackets in the first term on the right hand side is exactly the transpose ofS 0 from step 3: We conclude that the scheme in Algorithm 4.1 and the scheme in Algorithm 4.17 coincide in exact arithmetic, provided the matrix S 1 is invertible. However, the numerical behavior of the two schemes differs when S 1 is singular or close to singular. ForM n+1 close to singular, solving equation (25) might lead to numerical instabilities. This problem seems to be avoided in the projector-splitting scheme from [29,30], as no matrix inversion is involved. Such ill conditioning is however hidden in performing step 3 of Algorithm 4.17, since the QR or SVD decomposition can become unstable for ill-conditioned matrices (see [17, chap. 5]). In the case of a rank deficient basis {Ũ n+1 }, Algorithm 4.1 updates the stochastic basis by solving equation (25) in a least square sense while minimizing the norm Ỹ n+1 − Y n L 2 ρ . The previous subsection showed that such solution satisfies the discrete variational formulation which plays a crucial role in stability estimation (see Section 5.3). On the other hand, Algorithm 4.17 relies on the somehow arbitrary completion of the basis {U 1 } in the step 3. In presence of rank deficiency, the two algorithms can deliver different solutions (see Section 7.3 for a numerical comparison).

Stability estimates
The stability of the solution of problems similar to (6) are well analyzed (see e.g. [14]). A natural question is to what extent constraining the dynamics to the low rank manifold influences the stability properties. In Section 5.1, we will first recall some stability properties of the true solution u true of problem (6). Then, in Section 5.2 we will see that these properties hold for the continuous DLR solution as well. It turns out that our discretization schemes satisfy analogous stability properties, as we will see in Section 5.3. In particular, we will show that the implicit and semi-implicit version are unconditionally stable under some mild conditions on the size of the randomness in the operator. We will state two types of estimates: the first one holds for an operator L as described in Section 2 and a second one additionally assuming the operator L to be symmetric. Note that in the second case the bilinear coercive form ·, · L,ρ is a scalar product on L 2 ρ (Γ; V ). In the rest of this section we will assume that a solution of problem (6), a continuous DLR solution and a discrete DLR solution exist.

Stability of the continuous problem
We state here some standard stability estimates concerning the solution u true of problem (6).

Stability of the continuous DLR solution
Constraining the dynamics to the R-rank manifold does not destroy the stability properties from Proposition 5.1. Proof. Part 1: with 0, Y i L 2 ρ = 0, we can take u as a test function in the variational formulation (19). The rest of the proof follows the same steps as in the proof of Proposition 5.1.
. Asu ∈ V we can consideru as a test function in the variational formulation (19) and arrive at the sought result.
Part 3 and 4 are obtained analogously.

Stabilty of the discrete DLR solution
Now we proceed with showing stability properties of the fully discretized DLR system from Algorithm 4.1 for the three different operator evaluation terms corresponding to implicit Euler, explicit Euler and semi-implicit scheme. For each of them we will establish boundedness of norms and a decrease of norms for the case of zero forcing term f . The following simple lemma will be repeatedly used throughout.

Implicit Euler scheme
Applying an implicit operator evaluation, i.e. L(u n h,ρ , u n+1 h,ρ ) = L(u n+1 h,ρ ) results in a discretization scheme with the following stability properties.
for any time and space discretization parameters t, h > 0 with C L , C P > 0 the coercivity and continuous embedding constant defined in (4), (3), respectively.
In particular, for f = 0 and n = 0, . . . , N − 1 it holds: h,ρ L,ρ ≤ u n h,ρ L,ρ . Proof. Thanks to Theorem 4.9, we know that the discretized DLR system of equations with implicit operator evaluation can be written in a variational formulation as n = 0, . . . , N − 1.
1. Based on Lemma 4.6 we take v h = u n+1 h,ρ as a test function in the variational formulation (62). Using Lemma 5.3 results in Using the coercivity condition (4) and summing over n = 0, . . . , N − 1 gives us the sought result.

Explicit Euler scheme
Concerning the explicit Euler scheme (see subsection 4.1), which applies the time discretization L(u n h,ρ , u n+1 h,ρ ) = L(u n h,ρ ), the following stability result holds.
Theorem 5.5. Let {u n h,ρ } N n=0 be the discrete DLR solution as defined in Algorithm 4.1 with L(u n h,ρ , u n+1 h,ρ ) = L(u n h,ρ ). Then the following estimates hold: 2. If L is a symmetric operator we have Here C L , C B , C P > 0 are the coercivity, continuity and continuous embedding constants defined in (4), (5), (3), respectively and C I is the inverse inequality constant introduced in (22). For f = 0 and n = 0, . . . , N − 1 it holds: Proof. Thanks to the Theorem 4.9 we can rewrite the system of equations in the variational formulation 1. Based on Lemma 4.6 we take v h = u n+1 h,ρ as a test function in the variational formulation (65) and using Lemma 5.3 results in We further proceed by estimating where, in the third step, we used the inequality which holds based on assumption (22). Combining the terms, using the condition (63) and summing over n = 0, . . . , N − 1 finishes the proof.
2. Lemma 4.6 enables us to take u n+1 h,ρ − u n h,ρ as a test function in (65). This results in Using Lemma 5.3 we obtain where, in the second step, we used the assumption (22) 4. The proof of the forth property follows the same steps as the proof of part 2. Since there is no need to use the Young's inequality in (67), the condition on t/h 2p is weakened: As for the estimate in the · H,L 2 ρ -norm we can derive where in the last inequality we applied u n+1 h,ρ L,ρ ≤ u n h,ρ L,ρ for t h 2p ≤ 2 C 2 I C B .

Semi-implicit scheme
This subsection is dedicated to analyzing the semi-implicit scheme introduced in subsection 4.1 which applies the discretization L(u n h,ρ , u n+1 h,ρ ) = L det (u n+1 h,ρ ) + L stoch (u n h,ρ ). Apart from the inverse inequality (22) we will be using two additional inequalities. Let us assume there exists a constant C det > 0 such that This constant plays an important role in the stability estimation as it quantifies the extent to which the operator is evaluated implicitly. Its significance is summarized in Theorem 5.6. In addition we introduce a constant C stoch that bounds the stochasticity of the operator Theorem 5.6. Let {u n h,ρ } N n=0 be the discrete DLR solution as defined in Algorithm 4.1 with L(u n h,ρ , u n+1 h,ρ ) = L det (u n+1 h,ρ ) + L stoch (u n h,ρ ) with L det and L stoch satisfying (68) and (69), respectively. Then it holds 1.
2. If L is a symmetric operator we have Here C L , C B , C P , C I > 0 are the coercivity, continuity, continuous embedding and inverse inequality constants defined in (4), with t, h satisfying a weakened condition Proof. The variational formulation of the discrete DLR problem from Algorithm 4.1 reads in this case 1. We will consider v h = u n+1 h,ρ as a test function in (74) and we derive Combining the terms and summing over n = 0, . . . , N − 1 finishes the proof.

Explicit Euler scheme
Applying the explicit Euler scheme in the operator evaluation for a random heat equation, i.e. L(u n h,ρ , u n+1 h,ρ ) = −∇ · (a∇u n h,ρ ), results in the following system of equations The stability properties stated in Theorem 5.5 part 2 and 4 hold under the condition
Note that the condition (68) is automatically satisfied for a random heat equation, since we have and inf x∈D,ξ∈Γā a ≥ amin amax > 0.
The system of equations (31)- (33) can be rewritten as For a further specified diffusion coefficient we can state the following stability properties.
we have the stability properties (71) and (72) for any t, h.
Proof. The conditionā(x) ≥ a stoch (x, ξ) for every x ∈ D, ξ ∈ Γ implies i.e. C det ≥ inf x∈D,ξ∈Γā a ≥ 1 2 . Together with Theorem 5.6 we conclude the result. Proposition 6.1 tells us that applying a semi-implicit scheme to solve a heat equation with diffusion coefficient as described in (82) results in an unconditionally stable scheme. This result as well as some of the previous estimates will be numerically verified in the following section.

Numerical results
This section is dedicated to numerically study the stability estimates derived for a discrete DLR approximation in Section 5. In particular, we will be concerned with a random heat equation, as introduced in (78), with zero forcing term and diffusion coefficient of the form (82). We will look at the behavior of suitable norms of the solutions of the discretization schemes introduced in Section 4.1. We will as well look at a discretization scheme in which the projection is performed explicitly to see how important it is to project on the new computed basisŨ n+1 in (25). As a last result we provide a comparison with the projector-splitting scheme from [29].  with λ the Lebesgue measure restricted to the Borel σ-algebra B([−1, 1]). In this case the conditions (77), (29) and (68) are satisfied with a min > 0.04, C det > 1 2 . The initial condition is chosen as u 0 (x, ξ) = 10 sin(πx 1 ) sin(πx 2 ) + 2 sin(2πx 1 ) sin(2πx 2 )ξ 1 + 2 sin(4πx 1 ) sin(4πx 2 )ξ 2 + 2 sin(6πx 1 ) sin(6πx 2 )ξ 2 1 .
= 10 sin(πx 1 ) sin(πx 2 ) + 4 3 sin(6πx 1 ) sin(6πx 2 ) + 2 sin(2πx 1 ) sin(2πx 2 )ξ 1 The spatial discretization is performed by the finite element (FE) method with P 1 finite elements over a uniform mesh. The dimension of the corresponding FE space is determined by h-the element size. For this type of spatial discretization we have the inverse inequality (79): Concerning the stochastic discretization we will consider a tensor grid quadrature with Gauss-Legendre points for the case of a low-dimensional stochastic space M = 2 and a Monte-Carlo quadrature for the case M = 10. The time integration implements the explicit scheme and the semi-implicit scheme described in subsection 4.1. We will consider the forcing term f = 0, i.e. a dissipative problem and time T such that the energy norm ( · L,ρ ) of the solution attains a value smaller than 10 −10 . Our simulations were performed using the Fenics library [2].

Explicit scheme
Since f = 0, the result in Theorem 5.5 predicts a decay of the norm of the solution   Figure 1 shows the behavior of the energy norm ( · L,ρ ) and the L 2 norm ( · H,L 2 ρ ) in 3 different scenarios: in the first scenario we set h 1 = 0.142, t 1 = 0.0018, i.e. the condition t 1 /h 2 1 ≤ K is satisfied and observe that both the energy norm and the L 2 norm of the solution decrease in time (see Figure 1(a)); in the second scenario, we halved the element size h 2 = h 1 /2 and divided by 4 the time step t 2 = t 1 /4 so that the condition (85) is still satisfied. The norms again decreased in time (Figure 1(b)); in the third scenario we violated the condition (85) by setting h 3 = h 1 /2 and t 3 = t 1 /3. After a certain time the norms exploded ( Figure  1(c)).
To numerically demonstrate the sharpness of the condition (85), we ran the simulation with 72 different pairs of discretization parameters h, t. The results are shown in Figure 2, where we depict whether the energy norm at time T is bellow 10 −10 , in which case the norm was consistently decreasing; or more than 10 4 , in which case the solution blew up. We observe that a stable t has to be chosen to satisfy t ≤ Kh 2 , which confirms the sharpness of our theoretical derivations.

M = 10
In our second example we will consider a higher-dimensional problem: M = 10 for which we use a standard Monte-Carlo technique with 50 points. We observe a very similar behavior as in the small dimensional case. Figure 3 shows that satisfying the condition t 1 /h 2 1 ≤ K with K = 0.085 results in a stable scheme

Semi-implicit scheme
We proceed with the same test-case with M = 10, same spatial and stochastic discretization, i.e. Monte-Carlo method with 50 samples and employ a semiimplicit scheme in the operator evaluation. Since the diffusion coefficient considered is of the form (82) and f = 0, Theorem 5.6 predicts u n+1 h,ρ L,ρ ≤ u n h,ρ L,ρ ∀h, t, ∀n = 0, . . . , N − 1.
We set the spatial discretization h = 0.142 and vary the time step t. We observe a stable behavior no matter what t is used, which confirms the theoretical result (see Figure 4). We report that the results for M = 2 with 81 Gauss-Legendre collocation points exhibited a similar unconditionally-stable behavior.

Explicit projection
The following results give an insight into the importance of performing the projection in a 'Gauss-Seidel' way, i.e. projection on the stochastic basis is done explicitly, Y n kept from the previous time step, while the projection on the deterministic basis is done implicitly, i.e. we use the new computedŨ n+1 (see Algorithm 4.1 for more details). For comparison we consider a fully explicit projection, i.e. Y n as the stochastic basis and U n as the deterministic basis. We use a semi-implicit scheme to treat the operator evaluation term as described in subsection 4.1. As shown in Figure 5, in all 3 cases the solution reaches the zero steady state, however, not in a monotonous way.   Figure 5: Behavior of the energy norm ( · L,ρ ) for 3 different time steps when treating the projection in an explicit way (orange) and in a semi-implicit way (blue). We used the semi-implicit scheme for the operator evaluation term. We see that, as opposed to a semi-implicit projection, with an explicit projection we do not obtain an unconditional norm decrease

Conclusions
In this work we proposed and analyzed three types of discretization schemes, namely explicit, implicit and semi-implicit, to obtain a numerical solution of the DLR system of evolution equations for the deterministic and stochastic modes. Such discrete DLR solution was obtained by projecting the discretized dynamics on the tangent space of the low-rank manifold at an intermediate point. This point was built using the new-computed deterministic modes and old stochastic modes. We found this projection property to be useful when investigating stability of the DLR solution. The solution obtained by the implicit scheme remains unconditionally bounded by the data in suitable norms. Concerning the explicit and semi-implicit schemes, we derived stability conditions on the time step, independent of the smallest singular value, under which the solution remains bounded. Remarkably, applying the proposed semi-implicit scheme to a random heat equation with diffusion coefficient affine with respect to random variables results in a scheme unconditionally stable, with the same computational complexity as the explicit scheme. Our theoretical derivations are supported by numerical tests applied to a random heat equation with zero forcing term. In the semi-implicit case, we observed that the norm of the solution consistently decreases for every time-step considered. In the explicit case, our numerical results suggest that our theoretical stability condition on the time step is in fact sharp. Our future work includes investigating if the proposed approach can be extended to higher-order projector-splitting integrators, or used to show stability properties for other types of equations.