Dynamical low-rank 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 low-rank approximation has proven its efficiency include, e.g., kinetic theoryEiL18; einkemmer2019quasi; PeMF20; PeM20; einkemmer2020low; EiHY21; EiJ21; guo2021low as well as uncertainty quantification FeL18; MuN18; MuNV20; SaL09; kusch2021DLRUQ. In both fields the high-dimensional 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 projector-splitting 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 projector-splitting 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 anti-symmetry of the original problemCeL21. 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. Higher-order moments, however, are dampened heavily by the unconventional integrator and the projector-splitting 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 low-rank 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 low-rank approximation is applied. Second, the dynamical low-rank 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 projector-splitting integrator yields oscillatory first-order moments, while showing a satisfactory approximation for second-order moments compared to the unconventional integrator?
Should dynamical low-rank approximation be derived for the matrix ODE which results from a discretization of the original problem? Or should dynamical low-rank 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 non-linear DLRA evolution equations, despite being a tool for linear problems. This mainly stems from the fact that non-linearities 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 projector-splitting integrator applied to uncertain parabolic problemskazashi2020stability. 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 low-rank approximation. In contrast to previously derived schemes, the resulting discretization for the projector-splitting 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 low-rank 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 projector-splitting 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.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
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 stochastic-Galerkin approximations for linear problems with uncertainty gottlieb2008galerkin; gerster2019discretized. Let us discretize the above equation in space using a finite volume approximation with Lax-Friedrichs numerical flux. The spatial domain is decomposed into grid cells with equidistant spacing . A semi-discrete method for the solution , where then takes the form
The Lax-Friedrichs numerical flux with input reads
Writing the scheme without the definition of numerical fluxes gives
To simplify notation, we rewrite the time update in matrix notation. Let us define the tridiagonal matrices with non-zero entries in the off-diagonals
Then, when collecting the solution in , the scheme (2) becomes
Using a forward Euler time discretization with time step size and using , gives the fully discrete scheme
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
Here, denotes the imaginary unit. Collecting the basis functions in the matrices
The choice of our ansatz will simplify this scheme, since
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. 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
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 component-wise. Therefore, we have
Hence, for the Frobenius norm, we obtain
Due to Parseval’s identity, we have
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 low-rank 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
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
Choosing the discretization proposed in Section 2.1, the update of the full problem is composed of the two substeps
Written more compactly as a single update, the scheme becomes
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
With the CFL number , this gives
2.4 Dynamical low-rank approximation
In the following, we give a short overview on dynamical low-rank 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 low-rank approximation. The first one chooses a low-rank approximation on the matrix solution of (3) and the second one chooses a low-rank 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
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
where denotes the right-hand side of the semi-discrete scheme (3). We use to denote the tangent space of at . The stated problem can be reformulated (KochLubich07, Lemma 4.1) as
where denotes the orthogonal projection onto the tangent space, which is given by
The evolution equation (16) is then split by a Lie-Trotter splitting technique, yielding
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 projector-splitting integrator
-step: Update to and to via
Determine and with
by performing a QR decomposition.
-step: Update to via
and set .
-step: Update to and to via
Determine and with by performing a QR decomposition.
The time updated solution is then given by . For more details on the matrix projector-splitting 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
Determine with and store .
-step: Update to via
Determine with and store .
-step: Update to via
and set .
Note that these two integrators take the semi-discrete 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 stochastic-Galerkin) system or the scalar transport equation (or uncertain linear equation). Starting at a system of the form , the low-rank solution ansatz is
Note that we now have basis functions and . The corresponding split equations (17) become
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
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 projector-splitting integrator
3.1 Discrete dynamical low-rank approximation
In the following, we apply the projector-splitting 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
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
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 gives
If we define , this simplifies to
Now, since we know that the Fourier transform of the projector-splitting integrator takes the form (29), we can now investigate its stability properties.
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