Kinetic equations are a class of model equations used to describe the statistical behavior of a large number of particles that follow the same physics laws. They have been widely used in many aspects of physics and engineering: the classical Boltzmann equation describes the dynamics of rarefied gases, the radiative transfer equation characterizes the behavior of photons, the neutron transport equation describes the dynamics of neutrons in nuclear reactors, and the Vlasov–Poisson system has been used to describe the motion of plasmas.
These equations, despite there being significant differences for various kinds of particles, share a similar structure. In all situations the dynamics is described in phase space and the solution is thus a distribution function that counts the number of particles at a particular time , location , and velocity . One major reason for kinetic equations being challenging is that they are posed in a higher dimensional space; this is different from most other physical models that describe the corresponding dynamics only on the physical domain, e.g. for a fluid model the solution only depends on .
During the past decade, numerous methods have been proposed to numerically solve kinetic equations and many of them succeed in reducing computational cost without sacrificing accuracy. The main body of work concentrates on developing fast solvers for the collision operators and overcoming the stability requirement via relaxing the time-discretization [24, 22], through either finding fast methods to treat the transport term [17, 44, 5, 12], or implementing the resulting discrete system efficiently [18, 13]
. However, reducing the complexity due to the high dimensionality is largely left unaddressed. The difficulty here is rather clear. To numerically solve kinetic equations, one needs to sample a certain amount of discrete points in each dimension, resulting in a large number of degrees of freedom. The numerical cost, meanwhile, is typically determined mainly by the total number of degrees of freedom.
This viewpoint was challenged in  (based on earlier work on very high dimensional systems in quantum mechanics). There the authors abandon the traditional approach. They propose a dynamical low-rank approximation that relies more on the intrinsic dimensionality of the manifold the solution lives on. This method is designed to follow the flow of the dynamics and to identify the important features in the evolution. Since the equation is projected onto a solution manifold of lower dimensionality, only the core information is preserved and some redundant information is thrown away. The numerical cost thus depends mainly on the intrinsic dimensionality of the dynamics, rather than the total number of discrete points. The method was first utilized to deal with coupled ODE systems  and it has been demonstrated that, for a range of problems, it can preserve the accuracy of the solution.
The success in dealing with ODE systems inspired the generalization to PDEs. More recently, a class of efficient numerical schemes were proposed for kinetic equations. In the past few years, the method was applied to the Vlasov–Poisson system [28, 15], the Vlasov-Maxwell system  and the classical Boltzmann equation  in both collisionless and strong-collisional regimes, and it is observed that both the sophisticated Landau damping phenomenon for the VP system [45, 41], and the incompressible Navier-Stokes limit for the Boltzmann equation are captured rather accurately with low cost. Different equations may require slightly modified dynamical low-rank algorithms to suit the specific structures of the equation, but the general approach is quite similar: one looks for the main features in the evolution and follows the flow of the equation projected onto the solution manifold of a certain low rank.
However, despite the strong intuition and the many promising numerical experiments, a mathematical error analysis is largely absent from all of those works (especially in the PDE setting). There are two major difficulties here. First, we often do not have results to show that the true solution is indeed approximately of low-rank. This point is more subtle than one might think. In fact, a Fourier expansion is also a low-rank approximation. But, a Fourier expansion uses a fixed set of basis functions, leading to unnecessary high ranks. In addition, showing low-rank by performing a Fourier decomposition relies on assuming smoothness of the solution, which for kinetic equations is certainly problematic. It should be noted, however, that for elliptic equations some results can be obtained .
The second difficulty is that we have to show that the numerical method actually captures the low-rank structure of the solution. For the ODE case an analysis has been performed in  and related works, but this analysis only applies to the non-stiff case and is thus not applicable to PDEs. The only mathematical analysis of a dynamical low-rank scheme in the PDE setting we are aware of is found in . This paper considers a special situation where the stiffness originates only from the linear differential operator whose corresponding flow is exactly computable within the low-rank formulation. This separation suggests a splitting procedure that decouples the stiff and non-stiff dynamics, which permits the convergence of the resulting method. This is not the situation for kinetic equations, and thus the analysis does not apply.
The goal of the current paper is therefore to fill this void. In particular, we will consider a multi-scale radiative transfer equation, essentially a linearized Boltzmann equation. For this equation we will propose a dynamical low-rank algorithm and theoretically justify the validity and efficiency of this method in the strong-collisional regime (small Knudsen number). More specifically, we will show that in the small Knudsen number limit, the method is asymptotic preserving (that is, it correctly approximates the diffusion limit of the radiative transfer equation) and that it automatically captures the corresponding low-rank structure of the solution. It implies that the required degrees of freedom only scale as the number of grid points in physical space, in contrast to the full phase space as is required by other methods. This is, to our knowledge, the first numerical analysis result of this kind. Furthermore, our hope is that the present work will serve as a stepping stone for future studies for more complicated kinetic equations.
1.1. Multi-scale radiative transfer equation
We first give a quick overview of the equation. In dimension, the equation writes:
where is the scattering cross-section, , is called the Knudsen number, which represents the ratio between the mean free path and the typical domain length, and is the angular variable. The density is defined as follows
where is the area of , the unit sphere in dimensional space, and denotes the average with respect to . The equation is written in diffusion scaling, meaning the transport term and the collision term are multiplied by and , respectively.
Under certain regularity assumption on , in the limit of , the solution of (1) will asymptotically approach and thus loses its dependence. Then the solution of satisfies the following heat equation (with being a constant depending on ):
This is called the diffusion limit of the radiative transfer equation. This limit is particularly interesting in the present setting as the corresponding solution , indicating the solution is essentially of rank .
There are two major problems that pose difficulties for numerically solving equation (1). First, a small Knudsen number , the regime of interest in this paper, implies that both the transport term and the collision operator are stiff. A standard numerical scheme thus requires time steps that scale with and a very fine mesh. Second, equation (1) is posed in an dimensional phase space. Therefore, if the equation is discretized using points in each spatial direction and points in each velocity direction, we need to store at least floating point numbers. For this is a five-dimensional problem. The large increase in the degrees of freedom for such high-dimensional problems is usually referred to as the curse of dimensionality in the literature. Both of these issues contribute to the fact that solving equation (1) is extremely expensive from a computational point of view.
The second problem has been essentially left open. It is a rather widely-accepted fact that the numerical cost depends on the number of grid points (the degrees of freedom) used in the whole domain. However, there has been a large body of work addressing the first problem; namely, can we design a numerical solver that advances in time with a time-step size decoupled from the stiffness of the equation. A typical solution is to include some form of implicit solver in the numerical treatment in order to enlarge the stability region, and then devise a numerical method that efficiently solves the implicit part of the scheme. This property was later termed asymptotic-preserving (AP) in a ground-breaking paper , although some earlier schemes exists that were designed to satisfy this property . A systematic investigation of AP was then conducted in a series of works for many (mainly nonlinear) Boltzmann-like equations, see [8, 30, 3, 9, 31, 7]. Also see reviews [20, 23, 10].
1.2. Dynamical low rank approximation
Dynamical low rank approximation is a systematic approach to tackle the curse of dimensionality for time-dependent problems. Under the assumption that the solution in fact lives in a low-dimensional manifold, the method looks for the low-rank approximation to the solution at every time step and advances the evolution projected on a low-dimensional tangential space. In particular, for kinetic equations, it is rather straightforward to separate the physical space coordinate and the velocity space:
so that the low rank factors and depend only on one type of coordinate. Assuming , this approximation only requires degrees of freedom.
Historically, dynamical low-rank approximations have been considered extensively in quantum mechanics, in which the dimensionality is high. Finding a low rank approximation there makes the computation tractable, see [40, 39] and [33, 34, 4] for a mathematical treatment. The application of the method in a general setting is studied in [26, 27, 37, 1]
, where the authors study both the matrix case (the approach we will be using here), and a general tensor formats, as is required for very high dimensional problems from quantum mechanics. One significant disadvantage of these algorithm is that they are not robust with respect to over-approximation, i.e. choosing a rank larger than is necessary for the problem at hand. Initially this problem was solved by regularization, and then a major improvement was suggested in, in which the so-called projector splitting integrator
was introduced to make the dynamical low-rank approximation robust with respect to small singular values. This approach was later extended to various tensor formats [34, 35, 19, 38], and it is the approach to be utilized in the present paper.
The application of dynamical low-rank approximation for kinetic equations is relatively recent. Specifically, numerical methods have been proposed for the Vlasov–Poisson [15, 14], the Vlasov–Maxwell , and the Boltzmann equation . These are all nonlinear kinetic models and the numerical results are very promising. However, the currently available results all mainly focus on computational aspects and on conserving the physical structure of the underlying equation. A mathematical analysis is lacking.
To analyze the numerical error for a dynamical low-rank algorithm applied to kinetic equations, we have to answer two questions
Do the solutions have low-rank structures?
How can we capture the structure dynamically in numerics?
Intuitively, to answer the first question, we utilize the fluid limit obtained on the theoretical level going back to equations (1) and (3). Since one can show that as we get , the velocity direction completely degenerates, and the rank of the representation in equation (4) is simply . It is then reasonable to expect that for small the solution is only slightly different from its rank- approximation.
To answer the second question, however, requires us to design an appropriate algorithm. As mentioned above, we will use the projector-splitting approach. This leads to a set of three evolution equations for , , and , respectively. To complete one time step we have to advance all three equations. As will be shown in Section 3, however, the order in which the sub-flows in the splitting are solved is important for preserving the rank. In addition, we will demonstrate that the chosen time integrator plays a crucial role. One can, for example, show that the implicit Euler method, due to the lack of the symmetry, fails to satisfy the correct limit. On the other hand, the Crank–Nicolson method converges to the correct limit and is thus asymptotic preserving.
Before proceeding, let us put the present work in the proper context. There are two basic methods to compute a low-rank approximation of an evolutionary partial differential equations. The first option is to use a so calledprojection method (not to be confused with the projector-splitting integrator for a dynamical low-rank approximation). In this case we first discretize the partial differential equations. We then start with an initial value of fixed rank and compute one time step using the numerical integrator chosen. This, in general, takes us outside of the approximation space (i.e. the rank increases). The so obtained result is then projected back to a manifold of functions with fixed rank. For kinetic equations such an approach was suggested in . The disadvantage of this approach, however, is that an object that does not lie in the approximation space has to be constructed, which implies both a memory as well as a performance penalty. In addition, the derived algorithm is very closely tied to the specific space and time discretization that has been chosen. An alternative is to directly formulate the dynamics of the low-rank approximation by projecting the original equation onto the manifold of functions with a fixed rank. This can be done completely on the continuous level and results in a new set of partial differential equations formulated directly in terms the low-rank factors. As mentioned above, we will follow the latter approach which is usually referred to as a dynamical low-rank approximation.
The rest of the paper is organized as follows. In Section 2 we present the numerical method that is applied to the radiative transfer equation. In Section 3 we state and prove the main results. Some parts of the proof in Section 3 are rather tedious. We leave those to the appendix and keep only the core analysis in the main text. Numerical results are then presented in section 4.
2. Numerical scheme
In this section we follow the framework and notations in  and present a numerical method that uses a function of the form given in equation (4) as the approximation space. Our goal is to obtain a low rank approximation to the solution of (1). To define the rank, we first need to equip the spaces with proper measures. For that, we simply use the standard vector space in both the spatial and the velocity domain. That is, we have
where is the spatial measure in , and is a normalized measure in . Then, any function that can be expanded by a set of orthonormal basis functions in and , is called rank-. We collect all these functions together and denote the collection by .
Definition 1 (Rank- function in ).
The collection of all rank- functions is denoted by
where we call a function rank- if there is a set of orthonormal basis functions and a set of orthonormal basis functions such that
Here are orthonormal in physical space and are orthonormal in velocity space with appropriate measures and , namely
We note that in this definition, only the rank, , is fixed. The basis functions and can be arbitrary, as long as the orthogonality condition is satisfied. We further emphasize that is not a function space; it is easily seen that the summation of two rank- functions may not be rank-.
It is unlikely that the solution is of rank-, i.e. in , for all time. However, numerically one can argue that the solution is approximately of low rank. Thus for the numerical solution we seek a rank- approximation in at every time step. The algorithm is consequently looking for a trajectory on the manifold that resembles the evolution guided by the equation. In some sense, we need to project the equation in to the manifold and find the equation that governs the dynamics of this trajectory in . Let us denote by the analytic solution and by the numerical rank- approximation. Using the argument above, the governing equation for is
where stands for the projection of onto the tangential plane of at . This is to ensure that
where denotes the tangential plane of at . This condition is sufficient to guarantee that lies in the manifold for all time.
To explicitly express the tangential plane for the manifold and the projection operator, we first notice that at any function , the tangential plane is different and that the form highly depends on . Denote
then the tangential plane is given by
This set collects all functions whose infinitesimal, when added to , still yields a rank- function. In this definition, we notice that we are allowed to choose arbitrarily an matrix , a function list and a function list , as long as the gauge conditions, , are satisfied. We note that the gauge conditions are imposed to guarantee the uniqueness of the low-rank factors. Interested readers are referred to  for details.
With this definition, one has:
If , the orthogonality and gauge condition quickly allows us to write
if , its projection onto can also be easily formulated as follows
where the spatial and velocity projection are
The numerical method will then be developed upon this formulation. We discuss the semi-discrete (in time) and the fully-discrete schemes in details in the following subsections. This will give us a formulation for one time step. To devise an asymptotic preserving scheme for integrating from to the final time we need a thorough understanding of the error analysis. In fact, the implicit Euler and Crank-Nicolson behave slightly different. To optimally suit our purpose we run the Euler for one step before shifting to Crank-Nicolson scheme. We defer a detailed discussion to the end of Section 3.
We wrap up this introduction with a comment on rank- approximations. In fact, all numerical methods to solve PDEs are rank- approximations. Let us suppose we have a numerical method with grid points in both and . Then
where are basis functions we use to approximate the solution. For example, in a finite difference method we set as the Kronecker delta function or hat functions that peak at (similar argument holds for ) and is the numerical solution evaluated at and . To find a numerical solution on this mesh is equivalent to finding a rank- approximation with fixed and . Then plays the role of above. This is typically not formulated in this way as the rank is very large. When we consider a low-rank approximation, as in this paper, we will always assume that . To achieve this, the functions and need to evolve in time according to the dynamics of the equation.
2.1. Semi-discrete low rank splitting method
In this section, we develop a projector-splitting method to solve equation (7). In order to be concise, we simply denote the low rank solution by and we decompose the solution using its low-rank representation:
where and collects the basis functions
This formulation will also be quite useful later when we introduce a space discretization. We will also be using quantities and :
Computing (7) at discrete times amounts to finding the governing equations that provide the updates of
respectively, for all with .
To specify the initial data into the format, we project the initial condition onto
using the singular value decomposition (SVD):
For updating , the Lie-Trotter splitting is used. From time step to , we split the three operators on the right hand side of (7) into three sub-steps, namely
This splitting (12)-(14) takes place for each time step. That is, given the numerical solution at time we update the solution to by solving the three equations one after another. All equations are advanced for a full time step . Different from directly solving (7), each sub-step only changes one part of the decomposition. In particular, the first splitting step (12) preserves , the last (14) preserves , and the middle step (13) merely updates . This allows us to update the three components separately without disturbing others. Below we detail the evolution of each sub-step:
Starting with , in this sub-step, we run (12) for a full time step , and we denote the result by . Since the step preserves :
where , and and , both , are matrix versions of the differential operator and the scattering operator respectively:
We denote the solution to (16) by , and and
are obtained through the Gram-Schmidt process (QR factorization) which ensures the orthogonality of:
In this step, equation (13) is ran for a full time step with initial condition . We denote the result by . Since the step only changes , we immediately obtain
In matrix form this can be written as
where and , both , are the matrix versions of the multiplication operator associated to () and the density term, respectively:
We denote the computed update of (18) by .
In this step, (14) is ran for a full time step with initial value . We denote the result by . This step preserves . Thus,
This can be written in matrix form as follows
with (since )
Solving this equation we obtain , and the orthogonality of is ensured through the Gram-Schmidt process:
With these three steps completed, one finally arrives at numerical solution at time
2.2. Fully-discrete low rank splitting method
The beforementioned method will now be discretized. Due to the stiffness of the equations, an implicit scheme will be applied. We denote
the sets of discrete points in and . The discrete solution can then be represented as follows
with and denote the -th mode evaluated at the discrete points in and .
). The asymptotic analysis will be performed based on this fully-discrete formulation.
For updating (16), is preserved and thus . The direct application of the implicit Euler scheme gives
and if Crank–Nicolson is used, the scheme is written as
We have used
with . is an matrix with evaluations of at the grid points in assigned as diagonal entries and is the discrete approximation of . The specific form of depends on the spatial discretization. We have used the simple rectangle rule for the integration in . This determines the form of . Other numerical integral rules could also be applied and the specific form of will change accordingly. Obviously and are the discrete versions of (17).
Finally for updating (16), we note , , are preserved. Defining
and applying the implicit Euler method we obtain
where we used the same definition of . Crank–Nicolson method will not be used in the last step and thus we do not specify it.
We finalize the update for
3. Properties of the numerical scheme
We investigate the properties of the numerical method in this section. In particular we will discuss the computational complexity and prove the method preserves the asymptotic limit in the strong collisional regime.
3.1. Computational complexity
To analyze the computational complexity is rather straightforward. Denote the rank, and the number of grid points in physical space and velocity space respectively. The matrices in the updating formula, , , are computed by solving (22), (25), (29). The following cost incurs:
Preparation: Calculation of , needs floating point operations (flops).
Upadate : Calculation of needs flops.
Calculation of , need flops.
Solving needs flops.
Using QR decomposition to obtain , needs .
Update : Calculation of and needs
Calculation of , needs flops.
Solving needs flops.
Update : Calculation of needs needs flops.
Calculation of , needs flops.
Solving needs flops.
Using QR decomposition to obtain , needs .
Because , in conclusion, we need flops per time step, instead of the .
3.2. Intuition of the error analysis
Before describing and proving our results in detail, in this section we first give a relatively vague justification on why the method works. As described in the introduction, there are two points we need to make:
Why the true solution has an approximate low rank structure? This question is a rather fundamental, and is independent of the method chosen: the rank structure of the solution purely depends on the governing PDE we are studying here.
Why the method keeps track of the low rank structure? This question concerns the behavior of the specific method (dynamic low-rank approximation) we choose to use here.
The two questions will be addressed respectively in the following two subsections.
3.2.1. Foundation: the low rank structure in the solution
To justify the low rank structure of the linearized Boltzmann equation, we simply will cite the following result 
Theorem 1 (Theorem 2 of ).
Denote the solution to the linear Boltzmann equation
where and . Then in the zero limit of , the solution converges to the solution of the diffusion equation:
in the sense that
We note that periodic boundary conditions are imposed in the original theorem to avoid complications that may come from the boundary layers. Essentially this theorem states that in the zero limit of , loses its velocity dependence, and the dynamics will purely be reflected in the physical space. In some sense:
can be seen as the rank- approximation with the threshold set at any value bigger than .
The proof for the theorem follows the Hilbert expansion. Formally, writing , we plug it back in the original equation and get:
in which the last equation gives us the diffusion limit. As seen in the expansion, the essence of the proof mainly lies in showing that
and then tracing the dynamics of in the null space. Showing Theorem 1 rigorously then amounts to bounding uniformly in , and we omit it from here.
3.2.2. Capture the structure along the dynamics
It is not straightforward to show that the dynamical low-rank approximation method captures the rank structure of the solution. In fact, the direct calculation would suggest otherwise: recall the dynamic low-rank approximation solution governed by the following equation
and compare it with the original one (30), it is rather easy to see that one essentially need to show that
in some norm, where and are basis functions of . This does not seem to be an easy task. On the contrary, concerning that contains stiff terms such as , a brute-force calculation would suggest that the error is of order .
This straightforward approach, however, overlooks the information hidden in the equation. To show the method keeps small along the evolution, some delicacy from the equation needs to be employed. As will be presented in further details in Section 3.3, the proof largely relies on the separation of scales and one needs to perform the order by order matching of the scales to derive a clearer view of the solution structure. For that we perform the same analysis as done in (3.2.1), and asymptotically expand . Plugging it back into equation (32), one has, in the leading order:
Noting that from the definition of ,
which further suggests that the leading order of the numerical solution
This at least implies the reduced order equation (32) has its leading order lying in the correct space. Whether it preserves the correct equilibrium state is up to more delicate derivations in the higher order expansions, and thus is left to the following section, where we directly tackle the problem in the fully discrete setting, and show that the numerical scheme indeed captures the diffusion phenomenon.
3.3. Asymptotic preserving
Following the intuition from the previous section, we give the rigorous proof that shows that the method captures the diffusion limit in the limit.
It is a rather challenging task to design a numerical method for a stiff equation with the discretization being independent of the smallest scale. For an equation with small dependence, usually the solution shows variations at fine scale, and in order to preserve these variations, the discretization has to be small, leading to a large number of degrees of freedom, driving up the numerical cost. If a method has its discretization relaxed from requirements at the finest scale, while still preserves the asymptotic limit of the equation, we call it an asymptotic preserving (AP) method. It is an attractive property for a numerical method to have.
In this section, we will prove that in our method, by injecting the low-rank structure into numerical solutions, one can automatically capture the diffusion limit (3), with , , and all independent of , and thus is an AP method. Furthermore, can be as small as . We note that the analysis for Crank-Nicolson and implicit Euler shares some similarities but Crank-Nicolson enjoys a symmetry that allows us to pass to the asymptotic limit.
3.3.1. Asymptotic analysis of the implicit Euler method
To unify the notation, throughout this section, we denote a matrix of size with being the numerical solution at . Denote a column vector of length , then is a vector of length representing the discrete version of at time for spatial grid points . The hope is to show that density solves equation (3) in the limit .
Our first result concerns the behavior of the implicit Euler method in the limit .
We employ the implicit Euler method to compute