1 Introduction
Dynamical lowrank approximation (DLRA) KochLubich07
for parametrized partial differential equations has gained increasing attention in the last years. This stems mainly from its ability to mitigate the curse of dimensionality in terms of computational costs and memory requirements. Problems in which dynamical lowrank approximation has proven its efficiency include, e.g., kinetic theory
EiL18; einkemmer2019quasi; PeMF20; PeM20; einkemmer2020low; EiHY21; EiJ21; guo2021low as well as uncertainty quantification FeL18; MuN18; MuNV20; SaL09; kusch2021DLRUQ. In both fields the highdimensional phase space implies that obtaining numerical solutions is extremely expensive both in terms of memory and computational cost.Robust integrators for the DLRA evolution equations are the matrix projectorsplitting integrator, introduced in LubichOseledets, as well as the unconventional integrator, introduced in CeL21
. Both integrators are unaffected by the presence of small singular values. A main difference of the unconventional integrator is that the dynamics is only moving forward in time, whereas the projectorsplitting integrator includes a step, which moves backward. This property plays a key role in the stability for spatial discretizations, which we will see in this work. Furthermore, the unconventional integrator preserves symmetry or antisymmetry of the original problem
CeL21. On the other hand, the projector splitting integrator can be extended to second order, which has been widely used in the literature EiL18; einkemmer2019low; einkemmer2021asymptotic; einkemmer2019quasi.It has been shown in numerical experiments that the unconventional integrator yields smoother solution profiles for first order moments such as the scalar flux in radiation transport or the expected value in uncertainty quantification kusch2021DLRUQ. Higherorder moments, however, are dampened heavily by the unconventional integrator and the projectorsplitting integrator allows for a more adequate representation kusch2021DLRUQ. However, the reason for this behavior is not understood, which is the main motivation for this work.
In order to implement a dynamical lowrank integrator, the partial differential equation under consideration has to be discretized. There are two main approaches to determine an approximation to spatial derivatives in the equations of DLRA. First, the spatial discretization can be performed for the original equation, leading to a matrix differential equation to which the dynamical lowrank approximation is applied. Second, the dynamical lowrank approximation can be derived for the continuous problem and the evolution equations of DLRA can be discretized in space in a subsequent step (as has been suggested in EiL18). The first approach is extensively used. However, as we will show, it can suffer from instabilities. For the second approach, a set of differential equations is obtained that can be discretized by an appropriate method. This enables the implementation of a suitable stabilization for each individual substep of the two integrators. The construction of adequate stabilization strategies for each substep of the two integrators requires knowledge of dampening and amplification of the underlying dynamics, which we aim to establish in this work.
In this work, we answer the questions

Is there an analytic explanation why the projectorsplitting integrator yields oscillatory firstorder moments, while showing a satisfactory approximation for secondorder moments compared to the unconventional integrator?

Should dynamical lowrank approximation be derived for the matrix ODE which results from a discretization of the original problem? Or should dynamical lowrank approximation be performed for the continuous problem and the discretization be applied to the continuous DLRA evolution equations?

If the latter option is chosen: How should the time and space discretization for the different substeps be chosen to obtain a stable and accurate numerical method?
The tool that we use to answer these questions is a Fourier approach in the spirit of a von Neumann stability analysis. Remarkably, the Fourier analysis provides a deep understanding of the stability of the nonlinear DLRA evolution equations, despite being a tool for linear problems. This mainly stems from the fact that nonlinearities only arise in the basis functions, which, by Parselval’s identity, do not affect the
norm of the solution. The analysis recovers the behaviour seen in numerical experiments and enhances the understanding of dampening effects that are observed in the different DLRA approaches. To the best of our knowledge, a stability estimate of numerical schemes for DLRA evolution equations is only available in the case of the matrix projectorsplitting integrator applied to uncertain parabolic problems
kazashi2020stability. This analysis uses a discrete variational principle, which does not apply for hyperbolic and kinetic problems investigated in this work. The Fourier approach chosen in our work enables us to propose a stable discretization of the continuous projector splitting based dynamical lowrank approximation. In contrast to previously derived schemes, the resulting discretization for the projectorsplitting integrator is stable in the stability region of the full problem. Furthermore, we introduce a stable and efficient discretization of scattering for radiation transport. For this, we split scattering and streaming parts which is a common practice in radiation transport adams2002fast; hauck2013collision; crockatt2017arbitrary. By noting that the integrator for the scattering part only imposes dynamics in the step, we can omit the remainder, which reduces computational costs and provides a stable treatment of the scattering terms.This paper is structured as follows: After this introduction, we provide a general background to the used methods in Section 2 to give an overview on existing work and to fix notation. Here, we derive a spatial discretization of the original problem in Section 2.1, then we present a stability analysis of this discretization in Section 2.2 and include scattering into this scheme in Section 2.3. A short review of dynamical lowrank approximation is provided in Section 2.4, with a focus on the two robust integrators as well as their discrete and continuous formulations. Section 3 presents the stability analysis for the matrix projectorsplitting integrator. We start by pointing out potential stability issues for the discrete DLRA formulation in Section 3.1, and propose a stable discretization for the continuous formulation in Section 3.2. In Section 4, the proposed stability analysis is applied to the unconventional integrator which is shown to be stable even for the discrete DLRA formulation. An efficient and stable treatment of scattering terms that arise in kinetic transport is discussed in Section 5 and we provide numerical examples in Section 6.
2 Background
2.1 Discretization of the full problem
Parametric linear systems play an important role in various applications such as radiative transport or uncertainty quantification for material deformations. In the following, let us study a linear system of the form
(1) 
where we have and . Such systems can for example arise in the P or S approximations to radiative transfer case1967linear; lewis1984computational; adams2002fast as well as stochasticGalerkin approximations for linear problems with uncertainty gottlieb2008galerkin; gerster2019discretized. Let us discretize the above equation in space using a finite volume approximation with LaxFriedrichs numerical flux. The spatial domain is decomposed into grid cells with equidistant spacing . A semidiscrete method for the solution , where then takes the form
The LaxFriedrichs numerical flux with input reads
Writing the scheme without the definition of numerical fluxes gives
(2) 
To simplify notation, we rewrite the time update in matrix notation. Let us define the tridiagonal matrices with nonzero entries in the offdiagonals
Then, when collecting the solution in , the scheme (2) becomes
(3) 
Using a forward Euler time discretization with time step size and using , gives the fully discrete scheme
(4) 
2.2 stability analysis for the full problem
To recall certain details in the classical stability analysis and to fix notation, let us start by recalling the stability analysis for the full problem. Without loss of generality, we assume the spatial domain to be the interval . In this case, a discrete Fourier ansatz for the discretized solution takes the form
(5) 
Here, denotes the imaginary unit. Collecting the basis functions in the matrices
lets us write the Fourier ansatz (5) at time in matrix notation as . We use an upper case to indicate the adjoint matrix. The next step is to plug this wave ansatz into (4), which gives
The choice of our ansatz will simplify this scheme, since
(6) 
Here, the diagonal matrices have entries
Hence, the Fourier ansatz simplifies the scheme to
Now, with , we directly see that
We are interested in deriving an estimate for the Frobenius norm of (i.e. the
norm of the vector containing all degrees of freedom), which we denote by
. In the following we use to denote the spectral matrix norm and to denote the Euclidean norm for vectors.Collecting the Fourier coefficients in a vector gives the time update
Hence, when denoting the eigenvalue of as , the Euclidean norm gives for every
We thus have that
Hence, the eigenvalue which maximizes the amplification is , which denotes the biggest absolute eigenvalue of . Then, the amplification of a Fourier mode with wave number becomes
(7) 
Let us store the amplification factor in a diagonal matrix with
and collect the norm at wave number in a vector . Due to (7), the estimate holds componentwise. Therefore, we have
Hence, for the Frobenius norm, we obtain
Due to Parseval’s identity, we have
(8) 
When using the CFL number we obtain
To obtain stability we require an amplification factor which is smaller or equal to one, i.e., we must pick . For linear schemes, this stability (together with consistency) can be used to prove convergence. However, in this work, we focus on understanding dampening properties of the different integrators and leave the question of convergence for the (necessarily nonlinear) dynamical lowrank approximation to future research.
2.3 Stability for scattering terms
In the following, let us focus on the application of radiative transport. In this case, the original advection system (1) is augmented by scattering and absorption effects. This leads to the P equations, which read
(9) 
Here, denotes the Legendre polynomial of order and is a scattering matrix. The variable is the projected direction in which particles travel. To shorten notation, we define with entries . For isotropic scattering one for example has , i.e. scattering will not directly affect the scalar flux while dampening higher order moments. Let us investigate how the additional scattering affects stability. Commonly, scattering and streaming are treated separately through a splitting step, see e.g. adams2002fast. In this case, we can update the solution from time to time by
(10a)  
(10b) 
Choosing the discretization proposed in Section 2.1, the update of the full problem is composed of the two substeps
(11a)  
(11b) 
Written more compactly as a single update, the scheme becomes
(12) 
Let us again use a discrete Fourier ansatz . The next step is to plug this wave ansatz into (12), which gives
Now, with , we directly see that
Using Parseval’s identity yields
(13) 
with
With the CFL number , this gives
2.4 Dynamical lowrank approximation
In the following, we give a short overview on dynamical lowrank approximation KochLubich07 for problems of the form (1). The main idea of DLRA is to represent and evolve the solution on a manifold of rank functions. There are two approaches to derive the evolution equations of dynamical lowrank approximation. The first one chooses a lowrank approximation on the matrix solution of (3) and the second one chooses a lowrank approximation on the continuous level for the solution of the original problem (1), which is subsequently discretized.
Let us start by presenting DLRA for the discrete system (3). In this case, the solution is represented by
(14) 
where , and . The aim is to derive evolution equations for each of these factorization matrices. Let us denote the set of matrices that have the form (14) by . Then, we wish to find which fulfills
(15) 
where denotes the righthand side of the semidiscrete scheme (3). We use to denote the tangent space of at . The stated problem can be reformulated (KochLubich07, Lemma 4.1) as
(16) 
where denotes the orthogonal projection onto the tangent space, which is given by
The evolution equation (16) is then split by a LieTrotter splitting technique, yielding
(17a)  
(17b)  
(17c) 
This scheme can be used to update the solution from to . These split equations are reformulated to yield an efficient and robust integrator. Each substep in the above equations has a decomposition of the form (14). Defining the decompositions and gives the matrix projectorsplitting integrator

step: Update to via
(19) and set .

step: Update to and to via
(20) Determine and with by performing a QR decomposition.
The time updated solution is then given by . For more details on the matrix projectorsplitting integrator, we refer to LubichOseledets.
Recently, a further robust integrator, called the unconventional integrator, has been introduced in CeL21. This integrator works as follows:

step: Update to via
(21) Determine with and store .

step: Update to via
(22) Determine with and store .

step: Update to via
(23) and set .
Note that these two integrators take the semidiscrete matrix ODE system (3) as a starting point to derive DLRA evolution equations. I.e., the evolution equations are derived after performing the spatial discretization. Following EiL18, the DLRA evolution equations can also be derived for the continuous problem first, and the spatial discretization is performed on the DLRA equations second. Note that in our case, we are starting from a large system of partial differential equations (1), which can result from a discretization of the directional or uncertain domain of transport equations or linear equations with uncertainty. In our analysis, only the discretization of the spatial domain is important, which is why it does not matter whether the original problem is the P (as well as stochasticGalerkin) system or the scalar transport equation (or uncertain linear equation). Starting at a system of the form , the lowrank solution ansatz is
(24) 
Note that we now have basis functions and . The corresponding split equations (17) become
(25a)  
(25b)  
(25c) 
where we have and we choose to denote the inner product with respect to space. Furthermore, we will use the notation to indicate an integration over the spatial domain. Then, when collecting the spatial basis functions in the vector and storing the vectors as columns of the matrix , the corresponding , and equations read
(26a)  
(26b)  
(26c) 
Note that we use dots to indicate the time derivative for ordinary differential equations, whereas a partial time derivative is used for partial differential equations. The continuous formulation of the unconventional integrator takes the same
, and steps, but uses a different ordering. Note that the formulation (26) is continuous in space and requires a spatial discretization in order to evolve the system numerically in time. Compared to deriving the DLRA equations for the disrcete matrix ODE, this formulation provides more freedom in the choice of discretizations of each individual equation. At the same time, choosing such a discretization requires a profound understanding of the stability related to this set of equations.3 stability analysis for the matrix projectorsplitting integrator
3.1 Discrete dynamical lowrank approximation
In the following, we apply the projectorsplitting integrator to the matrix ordinary differential equation (3). Note that this corresponds to discretizing the full problem first and deriving the DLRA equations second. The , and steps from equations (18), (19) and (20) in combination with an explicit Euler time discretization then read
(27a)  
(27b)  
(27c) 
Here, we make use of the DLRA substeps , and . To underline similarities of the stability analysis to the full problem (cf. Section 2.2), let us go one step back to the corresponding split equations (17), which read
(28a)  
(28b)  
(28c) 
Omitting Roman indices, the solution of every substep in (28) is of the form , where . This is easily shown as every substep is of the form . We thus have . Therefore, one can choose , and . Then, the spatial discretization matrices and
can be Fourier transformed according to (
6), which givesIf we define , this simplifies to
(29a)  
(29b)  
(29c) 
Now, since we know that the Fourier transform of the projectorsplitting integrator takes the form (29), we can now investigate its stability properties.
Theorem 1.
Proof.
Let us pick a single mode solution with such that . In a more compact notation, we define the vector and with an arbitrary normalized vector , we have . Plugging this into the equations (28) yields for the first step
Here, we use that for our choice of the wave number we have and . Hence, the basis remains unchanged and only the coefficient changes its sign. Then for the second step, we have
Comments
There are no comments yet.