A low-rank Schwarz method for radiative transport equation with heterogeneous scattering coefficient

TL;DR

This paper introduces an accelerated Schwarz method leveraging random sampling to efficiently solve multiscale radiative transfer equations with heterogeneous scattering coefficients, demonstrating high accuracy and robustness.

Contribution

The work develops a novel low-rank Schwarz method using offline random sampling for local solution maps, improving efficiency in multiscale radiative transfer problems.

Findings

Method achieves high accuracy in numerical tests

Approach is robust across different heterogeneity levels

Significant reduction in computational cost

Abstract

Random sampling has been used to find low-rank structure and to build fast direct solvers for multiscale partial differential equations of various types. In this work, we design an accelerated Schwarz method for radiative transfer equations that makes use of approximate local solution maps constructed offline via a random sampling strategy. Numerical examples demonstrate the accuracy, robustness, and efficiency of the proposed approach.

Tables

Table 1. Table 1: Run time comparison between reduced Schwarz method with r = 2 , 3 , 4 , 5 , 6 𝑟 2 3 4 5 6 r=2,3,4,5,6 , an offline/online breakdown of Schwarz method and the vanilla Schwarz method.

	Running Time (s)
	$(\epsilon ,\delta )=(1/81,1/81)$		$(\epsilon ,\delta )=(1/81,1/9)$
	offline	online	offline	online
Low-Rank Schwarz $r=2$	81.01	0.0022	107.22	0.0024
Low-Rank Schwarz $r=3$	111.43	0.0023	148.09	0.0022
Low-Rank Schwarz $r=4$	139.97	0.0057	188.41	0.0024
Low-Rank Schwarz $r=5$	168.64	0.0032	227.99	0.0029
Low-Rank Schwarz $r=6$	197.49	0.0065	268.46	0.0162
Schwarz with full basis	535.26	0.0061	706.99	0.0148
Vanilla Schwarz	—	803.01	—	1027.40

Equations

$$v\cdot {\nabla }_{x}u(x,v)=\sigma (x)\mathcal{L}u(x,v),(x,v)\in \mathcal{D}:=\mathcal{K}\times \mathcal{V},$$

$$\epsilon v\cdot {\nabla }_{x}u(x,v)={\sigma }^{\delta }(x)\mathcal{L}u(x,v),$$

$$\Delta x\ll \min \{\epsilon ,\delta \}.$$

$$\epsilon v\cdot {\nabla }_{x}u(x,v)={\sigma }^{\delta }(x)\mathcal{L}u(x,v),$$

$$\mathcal{L}u(x,v)={\int }_{{\mathbb{S}}^{d-1}}u(x,v^{\prime })\mathrm{d}\mu (v^{\prime })-u(x,v),$$

$$k(x,v,v^{\prime })=\frac{1-g^2}{4\pi (1+g^2-2gv\cdot v^{\prime }{)}^{3/2}},$$

$$u(x,v)=u_0(x,v)+\epsilon u_1(x,v)+\cdots ,$$

$$\mathcal{O}(1):$$

$$\mathcal{O}(\epsilon ):$$

$$\mathcal{O}({\epsilon }^2):$$

$${\nabla }_x\cdot \left(\frac{1}{{\sigma }^{\delta }(x)}{\nabla }_x{u}_0(x)\right)=0,$$

$$\Vert u(x,v)-u_0(x){\Vert }_{L_2(\mathrm{d}x\mathrm{d}\mu (v))}=\mathcal{O}(\epsilon ).$$

$${\sigma }^{\delta }(x)=\sigma \left(x,\frac{x}{\delta }\right),$$

$$u^{\delta }(x,v)=u(x,y,v)=u\left(x,\frac{x}{\delta },v\right).$$

$$v\cdot {\nabla }_x{u}^{\delta }(x,y,v)+\frac{1}{\delta }v\cdot {\nabla }_y{u}^{\delta }(x,y,v)=\frac{\sigma (x,y)}{\epsilon }\mathcal{L}{u}^{\delta }.$$

$$u^{\delta }(x,y,v)=u_0(x,y,v)+\delta u_1(x,y,v)+\mathcal{O}({\delta }^2),$$

$$\mathcal{O}(1/\delta ):$$

$$\mathcal{O}(1):$$

$$i2\pi v\cdot\xi\hat{u}_{0}(x,\xi,v)=0\,,\quad\mbox{for all $\xi\in\mathbb{Z}^{d}$.}$$

$${\hat{u}}_0(x,\xi ,v)=0,\text{for almost all }(x,v).$$

$$u_0(x,y,v)=u_0(x,v),\text{for almost all }(x,v).$$

$$v\cdot {\nabla }_x{u}_0(x,v)=\frac{{\sigma }^{\ast }(x)}{\epsilon }\mathcal{L}{u}_0(x,v),\text{with}{\sigma }^{\ast }=\int \sigma (x,y)\mathrm{d}y.$$

$${\sigma }^{\delta }\stackrel{\delta \to 0}{\to }{\sigma }^{\ast },\text{in}{L}^{\infty }\text{weak-}\ast \text{topology},$$

$$\epsilon v\cdot {\nabla }_{x}u(x,v)={\sigma }^{\ast }(x)\mathcal{L}u(x,v).$$

$$\|\mathsf{A}-\mathsf{Q}\mathsf{Q}^{\top}\mathsf{A}\|_{2}\leq\left(1+17\sqrt{1+r/p}\right)\sigma_{r+1}+\frac{8\sqrt{k}}{p+1}\Bigl{(}\sum_{j>r}\sigma_{j}^{2}\Bigr{)}^{1/2}\,,$$

$$\epsilon v\cdot {\nabla }_{x}u(x,v)$$

$$u(x,v)$$

$${\Gamma }_{-}:=\{(x,v)\in \partial \mathcal{K}\times \mathcal{V}:-n_x\cdot v>0\}.$$

$$\mathcal{K}=\bigcup _{m=1}^M{\mathcal{K}}_m,$$

$$\mathcal{D}=\bigcup _{m=1}^M{\mathcal{D}}_m=\bigcup _{m=1}^M({\mathcal{K}}_m\times \mathcal{V}).$$

$${\mathcal{K}}_m^{\text{s}}:={\mathcal{K}}_m\setminus \left({\mathcal{K}}_{m-1}\cup {\mathcal{K}}_{m+1}\right),{\mathcal{D}}_m^{\text{s}}:={\mathcal{D}}_m\setminus \left({\mathcal{D}}_{m-1}\cup {\mathcal{D}}_{m+1}\right)={\mathcal{K}}_m^{\text{s}}\times \mathcal{V}.$$

$${\mathcal{E}}_{m,m-1}:={\mathcal{D}}_m\cap {\Gamma }_{m-1,-},{\mathcal{E}}_{m,m+1}:={\mathcal{D}}_m\cap {\Gamma }_{m+1,-},$$

Full text

1 Introduction

The radiative transfer equation (RTE) is a standard model that describes propagation of light through such turbid media as biological tissues or planetary atmospheres. The equation is used in situations in which energy is transported by light, as in the study of the greenhouse effect [4], optical tomography [28], and the radiation field for atmosphere-ocean system [38]. Light is injected from a source, and RTE models the absorption and scattering of the photons in the ambient material.

The model equation for the steady state is

[TABLE]

where $u(x,v)$ describes the light intensity at location $x$ oriented in velocity direction $v$. The left-hand side describes free propagation of the photons along direction $x$ with velocity $v$, while the right-hand side characterizes interaction between photons and media (via absorption and scattering). The media information is encoded in $\sigma(x)$, which is strictly positive for all $x$. The operator $\mathcal{L}$, typically an integral operator, characterizes how photons are scattered and change directions. We denote the physical domain by $\mathcal{K}$. Since photons always move with the same speed, the velocity term is determined purely by the direction, so that $v\in\mathcal{V}=\mathbb{S}^{d-1}$, the unit sphere in $d$ dimensions.

In the large space-regime, with scaling $x\to\frac{x}{\varepsilon}$ (where $\varepsilon$ is a small parameter discussed below), the equation (1) becomes

[TABLE]

where $\sigma^{\delta}(x)$ is the rescaled media function, with $\delta$ capturing the smallest scale of the variation in the media. This function is rough when $\delta\ll 1$. In the equation (2), $\varepsilon$ is the Knudsen number that represents the ratio of the mean free path to the typical domain length.

With appropriate boundary conditions, well-posedness of the equation is straightforward, and is independent of the scales (that is, the smallness of $\varepsilon$ or $\delta$) [2, 14]. In this paper, we tackle the numerical challenge of designing an efficient numerical solver for (2). We are especially interested in the regime of small $\varepsilon$ and small $\delta$, where classical numerical methods typically require high memory and computational cost, as we explain below.

1.1 Asymptotic preserving

Small parameters in PDEs can induce computational challenges. In the case described above, we have $\nabla_{x}u\sim{1}/{\min\{\varepsilon\,,\delta\}}$, so a classical numerical solver can be expected to attain good accuracy only when the mesh size in the discretization $\Delta x$ satisfies

[TABLE]

A grid in $d$ dimensions with this discretization parameter will have at least $N\gg\min(\varepsilon,\delta)^{-d}$ grid points, so the computation is prohibitive when $\varepsilon$ and $\delta$ are small.

A natural question is whether it is possible to design a numerical method for which the computational cost of obtaining a stable, accurate solution is independent of the parameters $\varepsilon$ and $\delta$, and whether the numerical solution can capture the right asymptotic limit of the solution as $\varepsilon$ and $\delta$ approach zero. If a numerical solver for a multiscale problem has its discretization independent of the smallest scale in the equation, but still preserves the asymptotic limits, then the solver is called asymptotic-preserving (AP). This term was coined in [26] for a class of kinetic equations, although some algorithms for simpler settings had been designed previously [30]. Extensive progress has been made during the past decade, with AP solvers being designed for the Bhatnagar-Gross-Krook equation (BGK, a special simplified version of the Boltzmann model that keeps the equilibrium), the Boltzmann equation, the Vlasov-Poisson-Boltzmann (VPB) equation, and many others [31, 13, 16, 15, 25]. See also the reviews [27].

A standard approach for designing AP solvers is based on analysis of the asymptotic limits. In some cases, asymptotic limits for the equations can be derived: the Euler limit for the Boltzmann, the coupled diffusion-Poisson system for the VPB system. One strategy for obtaining the AP property is to work with two sets of solvers, one for the original equation and one for the asymptotic limit, the latter being encoded in the former via a weight that can be tuned. In the limit as $\varepsilon\to 0$, this weight is adjusted so that the limiting equation solver dominates, driving the numerical solution to that of the asymptotic limiting system.

The analysis-based approach is straightforward and mathematically sound, and has made some previously impossible computations feasible. It depends, however, on analytical understanding of the asymptotic limit, which is not always straightforward. We are led to ask whether it is possible to design an AP solver that does not require detailed knowledge of the asymptotic limit. This paper addresses this question in the specific case of RTE. This equation is complicated in that different patterns of convergence of the pair $(\varepsilon,\delta)$ to $(0,0)$ lead to different limiting systems, not all of which are well understood. Can we design AP numerical solvers in the absence of this analytical understanding? We outline an answer to this question in the next section.

1.2 Random Sampling and PDE Compression

We design AP solvers without analytical knowledge by using compression techniques. Even when asymptotic limits of equations with small parameters are difficult to derive analytically, we can sometimes show the existence of such limits. In the discrete setting, a basic second order discretization scheme with $\Delta x=o(\sqrt{\tau}\epsilon,\sqrt{\tau}\delta)$ should suffice to attain the accuracy level $\tau$, which leads to $N_{\varepsilon}=\min\{{\varepsilon\sqrt{\tau}},\delta\sqrt{\tau}\}^{-d}$ grid points for $O(1)$ domain-size. When a limiting equation exists, this accuracy level may be attainable with as few as $N=\tau^{-d}$ grid points. Since the limiting equation is asymptotically close to the original equation, these $N$ degrees of freedom are asymptotically sufficient to represent the original PDE solution that naively would require $N_{\varepsilon}$ grid points to compute. This observation implies that the $N_{\varepsilon}$-dimensional solution space is compressible, and can be well approximated by a $N$-dimensional space when $\varepsilon$ and $\delta$ are small.

Knowing that the space is “compressible”, can one find the compressed space quickly? We answer this question affirmatively, in the case of the RTE, by making use of random sampling. Random sampling is not a new strategy. It has been used in data science to sample sparse vectors (as in compressed sensing [9]) and low-rank matrices [23], with the goal of reconstructing these objects from a relatively small number of samples. Generally, the number of samples is tied more closely to the intrinsic dimension of the object (for example, the number of non-zeros in a sparse vector) than to the dimension of the ambient space, which is typically much larger.

Applications of random sampling techniques to PDEs have been limited previously to the discrete algebraic systems obtained by discretizing the PDEs. A direct link to the original PDEs needs to be explored further. Some important questions have not been fully resolved, for example, whether the PDE solution space or the solution operator is “compressible”. The two views correspond to regarding a matrix as defining a column space or as a linear operator, respectively. Other questions involve the sense in which these objects are “compressible”, and whether the spectral norm used for matrices is the appropriate norm in the case of PDEs.

Previously, mostly in the context of elliptic PDEs, the low rank property of the solution space has been investigated and utilized in numerical solvers. For example, homogenization theory has been utilized [24, 37] for designing local basis functions for multiscale problems with structured media. For more general setting of $L_{\infty}$ media, the Kolmogorov $N$-width or the problem was studied in a pioneering paper [3], while the structure fo the Green’s functions were investigated in the framework of hierarchical matrix [6, 22]. Inspired by the studies, many algorithms have been proposed to utilize the rank (or decay) property, including [12, 35, 29]. Algorithms that specifically use PDE compression and random sampling ideas are developed in [8, 39, 36, 7, 17]. In the transport equation setting, [11] incorporated the random sampling technique within the discontinuous Galerkin framework for building local solution dictionaries. Corresponding to the first question asked above, in most of these papers, the authors regard the solution space to be “compressible” and the associated matrix is regarded as a column space. A more systematic investigation of PDE compression appears recently in our previous work [10], where compressed PDE solution spaces are related to low rank structure of the matrix formed by the Green’s functions. Such concepts in multi-scale PDE computation as asymptotic-preserving (see above) and numerical homogenization are unified under this framework.

1.3 Contribution

This paper follows the line of research started in [10]. For RTE (2), we know only that the equation has asymptotic limits with small parameters, but the actual forms of the limiting equations are unknown. We aim to design an accurate numerical scheme whose runtime is independent of the smallness of the coefficients in the equation.

We apply the Schwarz iteration under the domain decomposition framework. The domain is divided into overlapping subdomains (patches). The PDEs in these patches can be solved in parallel. Solution of the PDE on each patch with partial boundary conditions yields an output in the form of boundary conditions that are passed to neighboring patches. The PDE on each patch is solved again with the modified boundary conditions supplied by its neighbors, the whole process repeating until the solutions are consistent in the overlapping regions. The boundary-to-boundary map, in which the inputs are the partial boundary conditions on the patch PDEs and the output are the missing boundary conditions obtained by solving the PDEs, is a compressible map. We will develop an algorithm based on random sampling that computes an adequate approximation to this map quickly. The overall scheme is a composition of an offline component, in which low-rank approximations to the boundary-to-boundary maps are obtained using random sampling; and an online step, in which Schwarz iteration, accelerated by the low-rank boundary-to-boundary map, is executed until a solution consistent across the whole domain is found.

Our work contrasts with the approach in [10], where the local solution space is compressed in an offline step. In the online step, a solution for particular boundary conditions or source term is found as a linear combination of basis vectors for the compressed space, with the coefficients chosen to match the given conditions. The problem of finding these coefficients is typically overdetermined, the number of coefficients being fewer than the constraints arising from the boundary conditions or source term. Some accuracy is sacrificed, and the error is difficult to quantify. The current work compresses the boundary-to-boundary map, rather than the local solution space, in the offline stage, and uses the compressed map to update local boundary conditions in the online stage, until a preset error tolerance is achieved.

In this work, to demonstrate our numerical scheme and validate our theory, we consider a simpler setting of $1+1$ problem, that is, one spatial dimension and one velocity dimension. Our theory can be extended to higher dimensions in a conceptually straightforward way, but the implementation of the numerical scheme would become significantly more delicate in such cases. We do not pursue high-dimensional versions in this paper.

The remainder of the paper is organized as follows. We introduce the concept of “low-rankness” in the context of the RTE in Section 2. In Section 3, we review the Schwarz iteration under the domain decomposition framework, and present the new low-rank Schwarz iteration method based on random sampling. Numerical experience is described in Section 4.

2 Low-rankness of RTE in Various Regimes

As discussed above, current AP schemes rely heavily on good understanding of the analytical form of the asymptotic limits, although in some situations, this limiting form is hard to specify, even when we know that it exists. The radiative transfer equation with small Knudsen number $\varepsilon$ and small media oscillation period $\delta$ is a good example of the latter phenomenon. As $\varepsilon$ and $\delta$ converge to $(0,0)$ in different ways, the limiting equations are different, and only some of the limiting forms can be expressed explicitly. We show two different homogenization effects in the following two subsections, and unify them using the concept of the low-rankness in Section 2.3.

Consider the RTE in infinite domain (2), which we restate here:

[TABLE]

where $x\in\mathbb{R}^{d}$ and $v\in\mathcal{V}=\mathbb{S}^{d-1}$. We define the scattering operator $\mathcal{L}$ to have the following form:

[TABLE]

where $\mu(v)$ is the normalized measure on the velocity domain. The scattering coefficient $\sigma^{\delta}(x)>0$ encodes the media information, with $\delta$ denoting the smallest spatial scale. The operator has a nontrivial null space $\text{Null\,}\mathcal{L}$ which consists of functions that are constant in the velocity domain. We use this fact later to formally derive the diffusion limit of RTE. In this article, we choose the operator $\mathcal{L}$ to have this specific form, for simplicity. In practice, especially in applications to atmosphere science, the radiative transfer equation often takes this operator to be $\mathcal{L}u(x,v)=\int_{\mathbb{S}^{d-1}}k(x,v,v^{\prime})u(x,v^{\prime})\,\mathrm{d}\mu(v^{\prime})-\sigma(x,v)u(x,v)$ with the collision kernel

[TABLE]

where the constant $g\in[-1,1]$ determines the relative strength of the forward and backward scattering. This is the so-called Henyey-Greenstein model. The equation would have the same type of asymptotic limit (an elliptic equation) as long as $\sigma(x,v)=\int k(x,v,v^{\prime})\mathrm{d}\mu(v^{\prime})$ [2, 14].

We now consider different limits for different regimes of the parameters $(\varepsilon,\delta)$.

2.1 Diffusion Regime

In the diffusion regime, we have $\varepsilon\to 0$ while $\delta$ is fixed at a positive value. From (3), we see that $\mathcal{L}u(x,v)\sim 0$ in the leading order, meaning that $u(x,v)$ belongs to the null space of $\mathcal{L}$, and loses its velocity dependence. By matching orders in the classical asymptotic expansion

[TABLE]

we obtain

[TABLE]

By substituting the $\mathcal{O}(\varepsilon)$ equation into the $\mathcal{O}(\varepsilon^{2})$ equation, one obtains a diffusion equation, as follows.

Theorem 1** ([5]).**

In the zero limit of $\varepsilon$, the solution to (3) converges to the solution to the diffusion equation:

[TABLE]

in the sense that

[TABLE]

Remark 2**.**

*Note that we did not account for boundary conditions in deriving the limiting equation. In physical space, the derivation is valid when the boundary conditions are periodic. Otherwise, one has to be careful with the boundary influences and curvature effects. The diffusion limit still holds outside the boundary layers, but the convergence deteriorates when curvature corrections need to be taken into account. These results can be found in [21, 33] for the case when domain is convex. *

2.2 Homogenization Regime

When $\varepsilon$ is fixed at a positive value while $\delta\to 0$, homogenization limits are achieved; see [18]. We assume a two-scale media, having dependence on a fast variable $y=\frac{x}{\delta}$ and a slow variable $x$:

[TABLE]

where $\sigma(x,\cdot)$ is assumed to be periodic (with respect to the fast variable) for each $x$. Accordingly, we write the solution as

[TABLE]

In this notation, the operator $\nabla_{x}$ is replaced by $\nabla_{x}+\frac{1}{\delta}\nabla_{y}$ from chain rule. By substituting into the equation, we have

[TABLE]

By substituting the asymptotic expansion

[TABLE]

into the equation above, and matching terms, we obtain

[TABLE]

By applying the Fourier transform for the first equation (5a) with respect to the periodic variable $y$, we obtain

[TABLE]

We note that for almost all fixed $v\in\mathbb{R}^{d}$, the multiplier $i2\pi v\cdot\xi$ is non-vanishing for all $\xi\in\mathbb{Z}^{d}\backslash\{0\}$, because otherwise there exists some $\xi\in\mathbb{Z}^{d}\backslash\{0\}$ such that $i2\pi v\cdot\xi=0$ for a positive measure set of $v$, which is impossible. Therefore, by dividing the multiplier, we have for any $\xi\in\mathbb{Z}^{d}\backslash\{0\}$ that

[TABLE]

That is, all Fourier modes are vanishing except for the one with $\xi=0$. This implies that $u_{0}$ is independent of the periodic variable $y$, so we redefine the notation to omit this dependence:

[TABLE]

For the next order equation (5b), when we take the integral over $y$, the RHS vanishes due to periodicity, and we obtain the homogenized equation

[TABLE]

The derivation of homogenization limit presented above is validated in the following theorem. See [18, Theorem 3.1] for a rigorous proof for the time-dependent case.

Theorem 3**.**

Let $\sigma^{\delta}(x)$ be a bounded family of $L^{\infty}$ functions such that

[TABLE]

then the solution $u^{\delta}(x,v)$ to the RTE (3) (with $\varepsilon$ fixed) converges in $L^{\infty}$ weak-$\ast$ topology to $u(x,v)$, the solution to the following homogenized RTE:

[TABLE]

Because of the oscillations of scale $\delta$ in the media $\sigma^{\delta}(x)$, the solution $u^{\delta}$ is rough. However, such oscillations are homogenized in the $\delta\to 0$ limit, and the solution $u^{\delta}$ becomes close to the solution to (6), which has no oscillation.

In general, the limiting regime $\varepsilon\rightarrow 0,\delta\rightarrow 0$ can be taken through different routes. One may fix $\delta=1$ and send $\varepsilon$ to zero to reach the diffusion limit (4) and then send $\delta\rightarrow 0$, shown as solid arrow in Figure 1, Alternatively, one could fix $\varepsilon=1$ and send $\delta$ to zero to reach the homogenization limit (6), then send $\varepsilon\rightarrow 0$. This path is shown as the dashed arrow in Figure 1. Additionally, one could send both $\varepsilon$ and $\delta$ simultaneously to zero at different rates, shown as dotted arrows in Figure 1. All these routes, though considering the same regime $\varepsilon,\delta\rightarrow 0$, do not necessarily end up at the same limit. In fact, Goudon and Mellet [19, 20] showed that by following the route $\varepsilon=\delta\rightarrow 0$, RTE (2) ends up as an effective drift diffusion equation, while Abdallah, Puel and Vogelius [1] followed the route $\delta\gg\varepsilon\rightarrow 0$ to obtain an effective diffusion equation.

2.3 Low Rank of the PDE Solution Map

An AP scheme was proposed in [32] to deal with the regime $\delta\gg\varepsilon\rightarrow 0$, while numerical schemes for other regimes remain open. In practice, given a particular pair $(\varepsilon,\delta)$ that is close to zero, it is impossible to determine which limiting equation is the most appropriate one to use as an approximation to the solution of (2). The analysis-based approach of designing AP schemes is therefore not feasible. We seek to develop instead a universal numerical approach that is valid in different limiting regimes.

We start by considering the diagram in Figure 2.

Assume we are given an equation $\mathcal{R}^{\alpha}u^{\alpha}=0$ where $\alpha:=\min\{\varepsilon,\delta\}$ denotes the smallest parameters in the equation operator $\mathcal{R}^{\alpha}$, together with a boundary operator $\mathcal{B}$ such that $\mathcal{B}u^{\alpha}=f$ for some given boundary data $f$. The solution can be represented as a convolution of $f$ with all Green’s functions $\mathcal{G}^{\alpha}$, so it lies in the space spanned by $\mathcal{G}^{\alpha}$. To find an accurate numerical approximation to this solution, the operator $\mathcal{R}^{\alpha}$ is translated to $\mathsf{R}^{\alpha}$, a matrix with $N_{\alpha}\geq\frac{1}{\alpha}$ columns. Thus, the numerical solution $\mathsf{U}^{\alpha}$ is a vector of length $N_{\alpha}$.

On the other hand, assume the equation is “homogenizable” and there exists an asymptotic limit, an operator $\mathcal{R}^{\ast}$ so that the solution $u^{\ast}$ to equation $\mathcal{R}^{\ast}u^{\ast}=0$ with boundary condition $\mathcal{B}u^{\ast}=f$ is asymptotically close to $u^{\alpha}$. The computation of $u^{\ast}$ is expected to be significantly cheaper, since the limiting equation no longer has small parameters and is expected to be smooth. Thus, the numerical solution $\mathsf{U}^{\ast}$ requires merely $N=\mathcal{O}(1)$ degrees of freedom to represent $u^{\ast}$ accurately.

Since $u^{\ast}$ is close to $u^{\alpha}$ for regular $f$, the numerical solution spaces captured by the range of matrices $\mathsf{G}^{\ast}$ and $\mathsf{G}^{\alpha}$ are expected to be almost the same. On the other hand, although $\mathsf{G}^{\alpha}$ contains many more degrees of freedom (columns) than $\mathsf{G}^{\ast}$, the former is “compressible” and the latter is a good low-rank approximation to it. This argument can be made rigorous with the definition of “numerical rank of an operator,” a concept that is equivalent to the “Kolmogorov $n$-width”. More details can be found in [10].

2.4 Random Sampling for Low-Rank Structure

Knowing that an operator is approximately of low rank does not mean that it is easy to find a low-rank approximation quickly. In the linear algebra setting, finding the low-rank structure is equivalent to finding the singular vectors of a matrix that correspond to the largest singular values. For an $n\times m$ matrix, the singular value decomposition costs $\mathcal{O}(nm\min\{n,m\})$ operations, making the computation expensive for large matrices. When the approximate rank is known to be $r\ll n$, a randomized SVD (RSVD) solver based on a sketching procedure is available, whose cost depends on $r$. Properties of this approach are described in the following result.

Theorem 4** (Corollary 10.9 of [23]).**

Suppose that the matrix $\mathsf{A}\in\mathbb{R}^{n\times m}$ has singular values ordered as follows: $\sigma_{1}\geq\sigma_{2}\geq\cdots$. Assume that the target rank $r$ and oversampling parameter $p\geq 4$ are positive integers such that $k:=r+p\leq\min\{m,n\}$. Then with probability at least $1-6e^{-p}$, we have

[TABLE]

where $\mathsf{Q}\in\mathbb{R}^{m\times k}$ is a matrix with orthonormal columns whose column space matches that of $\mathsf{A}\Omega$, where $\Omega$ is a random matrix of dimension $m\times k$ with entries drawn from an independent identical distributed (i.i.d.) normal distribution.

This theorem suggests that if an operator has approximate low rank, random sampling can find its range accurately, with overwhelming probability. A simplified estimate shows that for oversampling parameter $p$ as small as $5$, with at least $99.8\%$ confidence, the error $\|\mathsf{A}-\mathsf{Q}\mathsf{Q}^{\top}\mathsf{A}\|_{2}$ can be controlled by $\left(1+11\sqrt{\min\{m,n\}(r+5)}\right)\sigma_{r+1}$. Algorithm 1, proposed in [23], finds the rank-$r$ approximation to $\mathsf{A}$, which we denote by $\mathsf{A}_{r}$.

There are two crucial features of the algorithm. First, the amount of computation depends crucially on the rank $r$, and is generally much less expensive than a full SVD. Second, it can be implemented without explicit knowledge of the matrix $\mathsf{A}$. Rather, we need only to be able to compute the products of $\mathsf{A}$ with the random matrix $\Omega$. These properties make the algorithm well suited for use in the numerical homogenization of PDEs.

RSVD is not the only algorithm that achieves the decomposition at the cost of $O(r)$ matrix-vector multiplications. Another approach is to explore a Krylov subspace of rank (slightly) greater than $r$ by initializing with a random vector and multiplying repeatedly by $\mathsf{A}$. The first $r$ singular values of the resulting matrix can be taken as approximating the leading $r$ singular values of $\mathsf{A}$. A variation of RSVD incorporates power iteration, which is similar to constructing a Krylov subspace; see [23, Sections 4.5 and 10.4]. We restrict here to the original RSVD algorithm for its simplicity and effectiveness in our numerical tests.

Translation of the randomization idea to the PDE setting is not straightforward. First, a PDE solution map is a continuous operator, not a matrix. Discretization of the space and the choice of norm is not always obvious. Redesigning the scheme to deal with an operator may also be difficult. While the adjoint of a matrix is easy to define, the adjoint of an operator may not be so easy to define. Second, a numerical homogenization scheme needs to find a solution quickly for an arbitrary boundary term or source, and in pursuit of that goal, the PDE solution map may not be the right operator to “compress”. In our approach, discussed in Section 3, we compress a different operator: the boundary-to-boundary map, Schwarz iteration scheme.

For the present, given an operator $\mathcal{A}:\mathcal{X}\to\mathcal{Y}$, we assume that the adjoint $\mathcal{A}^{\ast}$ is known. We also assume $\mathcal{X}$ and $\mathcal{Y}$ are finite dimensional and there is an inner product structure on $\mathcal{X}$ which allows us to efficiently draw random samples. We can then translate Algorithm 1 to the operator setting to find the corresponding Kolmogorov $r$-width operator $\mathcal{A}_{r}$ in Algorithm 2. (Note that upon fine discretization to meet the preset precision threshold, $\mathcal{X}$ can always be made finite dimensional.)

Remark 5**.**

Both Algorithm 1 and its operator counterpart Algorithm 2 require an input target rank $r$. Such parameter could be obtained from a priori estimate for elliptic type equations [3]. However, such guidance is unfortunately absent in the transport equation case. To address this issue, one may consider an adaptive randomized range finder (Algorithm 4.2 of [23]), in which a tolerance level is preset and the target rank is determined on the fly.

3 Low-Rank Schwarz Domain Decomposition Method

We consider the boundary value problem for RTE (3) in the following form:

[TABLE]

where the partial boundary $\Gamma_{-}$ is defined by

[TABLE]

Here, $n_{x}$ is the outer normal vector at location $x\in\partial\mathcal{K}$, and $v\cdot n_{x}$ is expected to be negative for all incoming velocities. Similarly, the outflow coordinates are collected in the complementary partial boundary $\Gamma_{+}:=\{(x,v)\in\partial\mathcal{K}\times\mathcal{V}:n_{x}\cdot v>0\}$. Problem (7) is well posed, as we show in the Appendix.

We start in Section 3.1 by introducing domain decomposition and the classical Schwarz method (Algorithm 3). Section 3.2 identifies the operator that needs to be compressed for efficient implementation of this method, while Section 3.3 derives the adjoint operator. These elements together make it possible to design the low-rank Schwarz method (Algorithm 5), which is presented in Section 3.4.

3.1 Schwarz Domain Decomposition Method

To solve (7), we first consider an overlapping domain decomposition of the physical space $\mathcal{K}$,

[TABLE]

where $\{\mathcal{K}_{m}\}_{m=1,2,\dotsc,M}$ forms an open cover of $\mathcal{K}$. We assume in the remainder of the discussion that the subdomains are ordered so that $\mathcal{K}_{m}$ can overlap only with $\mathcal{K}_{m-1}$ and $\mathcal{K}_{m+1}$. We decompose the domain $\mathscr{D}=\mathcal{K}\times\mathcal{V}$ accordingly as:

[TABLE]

We denote by $\Gamma_{m,\pm}$ the outflow and inflow boundaries for $\mathscr{D}_{m}$. We define those parts of the subdomains $\mathcal{K}_{m}$ and $\mathscr{D}_{m}$ that do not overlap with their neighbors as follows:

[TABLE]

Since the inflow boundary $\Gamma_{m\pm 1,-}$ of neighboring domain is partially inside $\mathscr{D}_{m}$, we further define

[TABLE]

so that the outflow boundary $\Gamma^{\text{s}}_{m,+}$ of $\mathscr{D}_{m}^{\text{s}}$ is the union of the domains above, that is,

[TABLE]

One can see that the restriction on $\Gamma^{\text{s}}_{m,+}$ of the local solution in the domain $\mathscr{D}_{m}$, would partially provide the inflow boundary condition for its neighboring domain $\mathscr{D}_{m\pm 1}$. This fact will be used later to update local solutions in each iteration of Schwarz method. Figure 3 illustrates our setup.

The Schwarz method is an iterative algorithm that updates solutions confined to different subdomains by exchanging information between iterations. In this setting, we update the values on $\mathcal{E}_{m,m\pm 1}$ using the newly computed solutions in each subdomain, repeating the process until the solution converges. To be more specific, we denote by $\phi^{k}_{m}$ the restriction of the solution at the $k$th step of the Schwarz process on patch $m$, confined to the partial boundary $\Gamma_{m,-}$, that is,

[TABLE]

The $k$th iteration of the Schwarz method can be expressed as a mapping from the $\phi^{k}$ to $\phi^{k+1}$, obtained by exchanging the boundary conditions between adjacent patches and solving the RTE. We denote this mapping by $\mathcal{P}$, as follows:

[TABLE]

We terminate at an iteration $k$ for which the difference between $\phi^{k}$ and $\phi^{k+1}$ falls below a given tolerance.

The evaluation of map $\mathcal{P}$ amounts to evaluation and assembly of the individual maps $\mathcal{P}_{m}(\phi_{m})$, $m=1,2,\dotsc,M$, and is defined by the following procedure.

• step 1

Define the following solution map

[TABLE]

by solving RTE on domain $\mathscr{D}_{m}$:

[TABLE]

Obtain the solution $u_{m}=\mathcal{S}_{m}(\phi_{m})$ in each subdomain $\mathscr{D}_{m}$ with boundary conditions $\phi_{m}$. Here $H_{A}(\mathscr{D}_{m})$ is a functional space where the trace of $u$ over the boundary $\Gamma_{m,\pm}$ is well defined (see the Appendix for details).

• step 2

Confine the solution $u_{m}$ on the boundaries:

[TABLE]

where $\phi^{+}_{m-1}$ and $\phi^{+}_{m+1}$ are the boundary values transmitted to the next iteration of the Schwarz procedure.

$\mathcal{P}_{m}$ is the boundary-to-boundary map that maps the boundary condition on $\Gamma_{m,-}$ to the boundary values on the adjacent subdomains ($\Gamma_{m,+}^{s}=\mathcal{E}_{m,m+1}\cup\mathcal{E}_{m,m-1}$):

[TABLE]

This map $\mathcal{P}_{m}$ is a well-defined operator, as we show in Theorem 36. It can be regarded as a composition of $\mathcal{S}_{m}$ and a trace operator, as follows:

[TABLE]

Note that $u_{m}|_{\Gamma^{\text{s}}_{m,+}}$ provides the boundary condition for the adjacent subdomains $\phi_{m\pm 1}$ in the next Schwarz iteration, seen as in equation (14). The full map $\mathcal{P}$ is obtained by collecting the action of $\mathcal{P}_{m}$ for all subdomains $m=1,2,\dotsc,M$.

As initial conditions for the Schwarz process, we set

[TABLE]

except at the physical boundary, where we impose given boundary conditions:

[TABLE]

When convergence to a given tolerance is achieved, the latest solutions may not perfectly match at the overlapping areas. To assemble the global solution, we define a suitable set of partition-of-unity functions $\{\eta_{m}(x)\}$ for the subdomains $\mathcal{K}_{m}$, whose properties are as follows:

[TABLE]

We construct the global solution by setting

[TABLE]

The method is summarized in Algorithm 3.

The Schwarz approach has several advantages. First, it is easy to implement in parallel, since the main computations (13) and (14) can be solved simultaneously for the subdomains $m=1,2,\dotsc,M$. In fact, one could even use different solvers in different subdomains, when appropriate (for example, when there is prior information about inhomogeneity of the medium). Second, computing solutions in each subdomain is significantly cheaper than for the full domain. It saves storage cost and computation time, especially when stoge and computation scale superlinearly with the size of the domain.

The disadvantage of the Schwarz approach is that it requires multiple iterations for convergence. Since $\mathcal{P}_{m}$ needs to be reevaluated at each iteration for each subdomain $\mathscr{D}_{m}$, and it calls for the computation of $\mathcal{S}_{m}$, finding the local solutions with the given boundary condition quickly is the key to the success of the entire algorithm. In the following sections, we identify the operator that can be efficiently compressed, aiming at improving the efficiency of evaluating $\mathcal{P}_{m}$, or $\mathcal{S}_{m}$.

3.2 Identifying the Operator to be Compressed

As discussed in Section 2.3, one should be able to reveal and exploit the low-rankness in homogenizable equations. In our setting, the local equation (13) has a homogenization limit when $\varepsilon$ and $\delta$ are small, so we expect the map from boundary conditions to local interior solutions (upon eliminating a boundary layer) to be of low rank. Indeed, we see this phenomenon in Figure 4, where we plot all normalized singular values of the discrete representation of $\mathcal{S}_{4}$ and $\mathcal{P}_{4}$. A solution with an inhomogeneous boundary condition can have strong boundary layer effect. These boundary layer effects are included in $\mathcal{S}_{m}$, an operator that maps the boundary condition to the solution in the entire region (including the boundary layer), destroying the desired low-rank structure. However, the operator $\mathcal{P}_{m}$ looks only at the solution confined to a small interior set $\Gamma^{\text{s}}_{m}$, and has a much faster decay in its singular values. This observation resonates with the argument in Remark 2: the homogenization limit concerns mainly the behavior of the solution in the interior of the (sub)domain, while the behavior of the solution in boundary layers is usually still far from “equilibrium”.

Thus $\mathcal{P}_{m}$ is a more suitable object for compression, so we seek a fast solver to approximate $\mathcal{P}_{m}(\phi_{m})$ for any input $\phi_{m}$, by making use of Randomized SVD (Algorithm 2). This method requires us to apply the operator $\mathcal{P}_{m}$ to random inputs, which amounts to finding the local solution in $\mathscr{D}_{m}$ with randomly constructed boundary conditions, and then confining it to $\Gamma^{\text{s}}_{m,+}$. Following (16), the evaluation $f=\mathcal{P}_{m}(r)$ is defined as

[TABLE]

Finding $\mathcal{P}^{\ast}_{m}$, the adjoint of $\mathcal{P}_{m}$, is much more complicated, as we show in the following theorem, whose proof appears in the Appendix.

Theorem 6**.**

The adjoint operator $\mathcal{P}^{\ast}_{m}$ is defined by

[TABLE]

where $h$, supported on $\mathscr{D}_{m}\backslash\mathscr{D}_{m}^{s}$, satisfies:

[TABLE]

in which $g$ is the solution to:

[TABLE]

The operators $\mathcal{P}^{\ast}_{m}$ and $\mathcal{P}_{m}$ are adjoint in the sense that:

[TABLE]

where $\langle\cdot,\cdot\rangle_{\Gamma_{m,+}^{s}}$ and $\langle\cdot,\cdot\rangle_{\Gamma_{m,-}}$ are weighted-$L^{2}$ inner products on $L^{2}(\Gamma^{\text{s}}_{m,+};|n\cdot v|)$ and $L^{2}(\Gamma_{m,-};|n\cdot v|)$, respectively, defined by

[TABLE]

The computation involved in finding the adjoint operator is complicated. It requires the computation of two adjoint RTEs over $\mathscr{D}_{m}^{\text{s}}$ and $\mathscr{D}_{m}\backslash\mathscr{D}_{m}^{\text{s}}$, respectively, that are coupled in a nontrivial fashion through the boundary conditions, as seen in (21) and (22). We further note that the measure is not the standard Lebesgue measure, but rather is weighted by $|n\cdot v|$.

3.3 Design of the Adjoint Map for $\mathcal{S}^{\text{s}}_{m}$

Although the operator $\mathcal{P}_{m}$ is of approximate low-rank, its adjoint $\mathcal{P}_{m}^{\ast}$, which is needed to compute the low-rank approximation, is complicated. The operator $\mathcal{S}_{m}$, on the other hand, is not compressible, but its adjoint is relatively easy to find. In this section, we show that we can approximate $\mathcal{P}_{m}$ by an approximately low-rank operator based on $\mathcal{S}_{m}$ whose adjoint is easy to find.

Since $\mathcal{S}_{m}$ has slow singular decay mainly because it contains too much information from the boundary layer, we consider a restriction of this operator from $\mathscr{D}_{m}$ to $\mathscr{D}_{m}^{\text{s}}$, which we call $\mathcal{S}^{\text{s}}_{m}$. The restriction to $\mathscr{D}_{m}^{\text{s}}$ eliminates most of the effects of the boundary layer. This operator is defined as follows.

[TABLE]

where $u_{m}^{\text{s}}=u_{m}|_{\mathscr{D}^{\text{s}}_{m}}$ and $u_{m}=\mathcal{S}_{m}\phi$. The advantages of using this operator are threefold.

$\mathcal{P}_{m}$ can be defined easily in terms of $\mathcal{S}^{\text{s}}_{m}$. Nothing is lost by comparison with (16); we have

[TABLE]

and $u^{\text{s}}_{m}|_{\Gamma^{\text{s}}_{m,+}}$ once again serves as the new boundary condition $\phi_{m\pm 1}$, as in equation (14). Note that the trace in (25) is well defined, as $\mathcal{S}_{m}$ maps boundary conditions to $H_{A}(\mathscr{D}_{m})$, so the image of its restriction to $\mathscr{D}_{m}^{\text{s}}$ has a trace on the boundary $\Gamma_{m,+}^{\text{s}}$ of $\mathscr{D}_{m}^{\text{s}}$.

• 2.

Because effects from boundary layers are excluded in $\mathcal{S}_{m}^{\text{s}}$, it can be expected to have approximate low rank. Figure 5 shows that the decay rate of $\mathcal{S}^{\text{s}}_{m}$ (upon discretization) is almost the same as for $\mathcal{P}_{m}$.

The adjoint is easy to compute, as we show next, in Theorem 7.

Theorem 7**.**

The adjoint of $\mathcal{S}_{m}^{\text{s}}$ is defined as follows:

[TABLE]

where $h$ solves the adjoint RTE over $\mathscr{D}_{m}$, which is

[TABLE]

and the source $\tilde{g}$ is the trivial extension of $g$ over $\mathscr{D}_{m}$, that is,

[TABLE]

Proof.

We need to show that

[TABLE]

for all $g$ and $\phi$. Denoting by $u$ the solution to (13) with boundary condition $\phi$, then the definition of $\mathcal{S}^{\text{s}}_{m}$ implies that the left hand side of this expression is

[TABLE]

Denoting by $h$ the solution to (27) with source term $g$, we have

[TABLE]

yielding the desired result. ∎

By comparing Theorem 6 with Theorem 7, we see immediately that computing the adjoint operator $\left(\mathcal{S}_{m}^{\text{s}}\right)^{\ast}$ is significantly easier than computing $\mathcal{P}^{\ast}_{m}$.

3.4 Low-Rank Schwarz Iteration Method

We can use $\left(\mathcal{S}_{m}^{\text{s}}\right)^{\ast}$ to implement the RSVD method to find the low-rank approximation to the operator $\mathcal{S}_{m}^{\text{s}}$. Given target rank $r$, and denoting the reduced operator by $\mathcal{S}_{m,r}^{\text{s}}$, we look for functions $\mu_{i}$ and $\nu_{i}$, and nonnegative scalars $\sigma_{i}$, such that:

[TABLE]

where $\mu_{i}(x_{1},v_{1})$ and $\nu_{i}(x_{2},v_{2})$ are obtained in Algorithm 4.

We note that in this algorithm, $k$ and $r$ could be hard to choose ahead of time. Numerically we can choose it “on-the-fly”. This means we simply set an accuracy threshold and stop the process once the newly computed $u_{m,i}|_{\mathscr{D}_{m}^{s}}$ falls almost in the previous generated space within the preset error tolerance. In terms of the low-rank operator $\mathcal{S}^{\text{s}}_{m,r}$, the procedure (25) is reduced further to

[TABLE]

and we once again use $u^{\text{s}}_{m}|_{\Gamma^{\text{s}}_{m,+}}$ to obtain the solutions $\phi_{m\pm 1}$ to be used at the next time step, as in equation (14).

The procedure we have just outlined provides a much cheaper way to evaluate $\mathcal{P}$ in (16) for the following reasons.

$\mathcal{S}_{m}^{\text{s}}$ maps the boundary condition to the interior of the subdomain, and it is cheaper to evaluate than $\mathcal{S}_{m}$, whose range has a bigger support. 2. 2.

The format in (28) guides the evaluation of $\mathcal{S}_{m,r}^{\text{s}}(\phi_{m})$; we have

[TABLE]

This evaluation requrires $\mathcal{O}(k|\Gamma_{m,-}|)$ operations, where $|\Gamma_{m,-}|$ is the cardinality (the number of grid points) in $\Gamma_{m,-}$.

Algorithm 5 summarizes the complete approach using the reduced solution map. The method is divided into offline and online stages. The reduced operators $\mathcal{S}^{\text{s}}_{m,r}$ are found in the offline stage, then called repeatedly in the online stage, during the Schwarz iteration procedure.

4 Numerical examples

In this section, we present numerical examples to validate the accuracy and efficiency of our methods. We consider boundary value problem (7) with domain $\mathscr{D}=\mathcal{K}\times\mathcal{V}=(0,1)\times(-1,1)$ and highly oscillatory scattering coefficient $\sigma^{\delta}(x)$ defined by

[TABLE]

where $\delta$ represents the period of oscillation in the spatial space. See Figure 6 for a graph of $\sigma^{\delta}(x)$ with $\delta=1/81$.

The space domain $\mathcal{K}$ is divided into $M=10$ different local subdomains $\mathcal{K}_{m},m=1,\ldots,M$, as follows:

[TABLE]

so that each subdomain $\mathcal{K}_{m}$ overlaps with its neighboring subdomains $\mathcal{K}_{m-1}$ and $\mathcal{K}_{m+1}$ (except for the subdomains $\mathcal{K}_{1}$ and $\mathcal{K}_{M}$ at the two ends of the domain, which overlap with just one neighbor each). The subdomains $\mathcal{K}_{m}^{\text{s}}$ defined in (10) are

[TABLE]

We thus have $\mathcal{E}_{m,m\pm 1}=\{\frac{m}{M}\}\times\mathcal{V}\,,\quad m=1,\ldots,M-1$, while $\mathscr{D}_{m}^{\text{s}}:=\mathcal{K}_{m}^{\text{s}}\times\mathcal{V}$ satisfies

[TABLE]

The spatial domain $\mathcal{K}$ is discretized with a extreme fine mesh of size $\Delta x=1/360$, while the velocity domain $\mathcal{V}$ is discretized with a mesh of size $\Delta v=2/40=0.05$. The fine-mesh discretization is determined by $\delta$ and $\varepsilon$.

4.1 Local Tests

We first show that the singular values of the local solution map $\mathcal{S}_{m}^{\text{s}}$ indeed decay rapidly for small Knudsen number $\varepsilon$ and small $\delta$. Figure 7 plots the singular values of $\mathcal{S}_{4}^{\text{s}}$ (relative to the largest singular value) for various values of $(\varepsilon,\delta)$. In the case of large values $\varepsilon=\delta=1$, the singular values decay slowly and low-rank structure is not present. By contrast, in the small-value regimes $(\varepsilon,\delta)=(1/81,1/81)$ and $(\varepsilon,\delta)=(1/81,1/9)$, low-rank structure is evident. As a consequence, only half or even a quarter of basis functions are needed to achieve high accuracy in approximating the local solution map $\mathcal{S}_{m}^{\text{s}}$.

4.2 Global Tests

We consider solving RTE (7) with scattering parameter (31) and the following inflow boundary conditions over $\Gamma_{-}$:

[TABLE]

We approximate $\mathcal{S}_{m}^{\text{s}}$ by low-rank operators $\mathcal{S}_{m,r}^{\text{s}}$, according to Algorithm 4, with $r=2,3,4,5,6$. We then use these low-rank approximations in the reduced Schwarz method, Algorithm 5. The various approximating solutions are then compared to the reference solution, also obtained by the Schwarz method, in terms of accuracy and speed, for different values of the parameter pair $(\varepsilon,\delta)$. We also document the global error as a function of the number of iterations.

Accuracy of approximating solution

Figures 8, 9, and 10 show the reference solution over domain $\mathscr{D}$ and compare with approximate solutions for $r=2$ and $r=6$, for parameter pair settings $(\varepsilon,\delta)=(1/81,1/81)$, $(1/81,1/9)$, and $(1,1)$, respectively. In Figure 8, for $(\varepsilon,\delta)=(1/81,1/81)$, the approximate solutions are very close to the reference solution. Figure 9, with $(\varepsilon,\delta)=(1/81,1/9)$, shows poor approximation for $r=2$ but good approximation for $r=6$. For the large-value case $(\varepsilon,\delta)=(1,1)$, shown in Figure 10, both approximations are poor, due to the lack of low-rank structure in $\mathcal{S}_{m}^{\text{s}}$.

Figure 11 shows the relative difference between approximate and reference solutions, plotted as a function of $r$, for the three settings of $(\varepsilon,\delta)$ considered here. The difference is defined by the formula

[TABLE]

where $u_{\text{approx}}$ is the approximate solution in question and $u_{\text{ref}}$ is the numerical reference solution computed with fine mesh. We evaluate the difference using $l_{2}$ norm of the two vectors. We see that the quality of the approximate solution aligns with the local singular value decay shown in Figure 7. For large value case $(\varepsilon,\delta)=(1,1)$, there is no decay in relative errors as $r$ increases. For $(\varepsilon,\delta)=(1/81,1/81)$, the relative error is below $10\%$ for $r=3$ and decreases as $r$ increases. For $(\varepsilon,\delta)=(1/81,1/9)$, the relative error decreases rapidly with $r$.

Efficiency of approximating solution

To demonstrate the efficiency of our method, we compare our reduced Schwarz method and the “vanilla” Schwarz method (which does not use low-rank approximations) in terms of accuracy and running time. In particular, we run the reduced Schwarz method for $T=50$ iterations, and compare the number of iterations needed for vanilla Schwarz method to achieve the same accuracy. Figure 12 plots relative error as a function of iteration number $t$ of reduced Schwarz method and the vanilla Schwarz, for parameter pairs $(\varepsilon,\delta)=(1/81,1/81)$ and $(\varepsilon,\delta)=(1/81,1/9)$ and rank $r=6$. The convergence speed of the vanilla and reduced versions of the Schwarz methods are quite similar. However, the reduced Schwarz iteration is significantly cheaper due to the use of the low rank structure. We document the run time on a standard PC in Table 1. We note that the vanilla Schwarz iteration does not have the offline step, and the online stage amounts to computing the equation for the given boundary condition on each subdomain, and is very expensive. If we need to solve the RTE for multiple boundary conditions, the computational saving would be quite significant. We also report on results obtained with the Schwarz procedure in which the full basis is prepared offline (denoted in the table as “Schwarz with full basis”).

5 Conclusion

Random sampling is a popular technique in data science to find low-rank structure of a matrix. In particular, the randomized SVD method has proved to be an efficient and accurate technique for finding the dominant singular values and singular vectors of a matrix at reasonable cost.

Partial differential equations with multiscale structures can usuallly be described by an effective equation without fine scale details. In some sense, this phenomenon means the operator and the solution space are of low rank. We exploit this property to design efficient numerical schemes, using the radiative transfer equation as an example. Specifically, we utilize the Schwarz iteration in the domain-decomposition framework, with a boundary-to-boundary map that communicates between subdomains. We use Randomized SVD to approximate this boundary-to-boundary map, exploiting the low-rank structure of the map. This computation is performed offline; the Schwarz iteration that makes use of these approximate maps is performed online. Numerical examples confirm the effectiveness and computational efficiency of our approach.

Several aspects and extensions of our approach remain to be investigated. The biggest obstacle in making our approach fully rigorous is the lack of theoretical guarantees on the decay of singular values of the PDE operator, in either the $\epsilon\to 0$ or the $\delta\to 0$ regimes. The lack of such guarantees makes numerical analysis hard to perform. Such results might point the way to better designs of the domain partition (especially the size of the buffer zones) and better choices of the accuracy threshold that defines the target rank $r$.

Extensions to spatial domains of dimension greater than $1$, especially in the choice of domain partitioning and computation of the boundary-to-boundary maps, remain a significant computational challenge that could be addressed in future work.

Acknowledgment

The authors thank the two anonymous referees for insightful suggestions, which lead to great improvement of the paper.

Appendix A Well-posedness Theory of RTE

We show the well-posedness theory of RTE with fixed parameters $\varepsilon,\delta>0$ in this appendix. Most results are cited from  [2]. For simplicity, we consider the following RTE with inflow boundary condition:

[TABLE]

First, we introduce a functional space $H^{1}_{2}(\mathscr{D})$ with norm defined as follows:

[TABLE]

To find a suitable functional space for solutions of RTE, we modify $H_{2}^{1}(\mathscr{D})$ and define a Hilbert space $H_{A}(\mathscr{D})$ with the following scalar product and norm:

[TABLE]

$H_{A}(\mathscr{D})$ is obviously a subspace of $H_{2}^{1}(\mathscr{D})$. The following theorem shows well-posedness of RTE over $H_{A}(\mathscr{D})$:

Theorem 8** (Theorem 3.7 in [2]).**

Given an inflow boundary condition $\phi\in L^{2}(\Gamma_{-};|n\cdot v|)$, then there exists a unique solution $u\in H_{A}$ to RTE such that

[TABLE]

Further, the trace operator $\gamma_{\pm}:u\in H_{2}^{1}(\mathscr{D})\mapsto u|_{\Gamma_{\pm}}\in L^{2}(\Gamma_{\pm};|n\cdot v|)$ is also well-defined by the following theorem:

Theorem 9** (Theorem 2.8 in [2]).**

If $u\in H_{2}^{1}(\mathscr{D})$, then $u$ has a trace over $\Gamma_{\pm}$ belonging to $L^{2}(\Gamma_{\pm};|n\cdot v|)$. In addition, we have

[TABLE]

Moreover, for RTE with a source term:

[TABLE]

the following theorem holds.

Theorem 10** (Theorem 3.13 in [2]).**

If $f$ belongs to $H_{A}^{-1}(\mathscr{D})$, the dual space $H_{A}$ under $L^{2}$ pairing over $\mathscr{D}$, $\phi\in L^{2}(\Gamma_{-};|n\cdot v|)$ then the above equation admits an unique solution $u\in L^{2}(\mathscr{D})$ such that

[TABLE]

Further, if $f\in L^{2}(\mathscr{D})$, then the solution $u\in H_{A}(\mathscr{D})$.

Remark 11**.**

The theorems above also hold for adjoint RTE with or without a source term. Details can be found in [2].

Proof of Theorem 6

Proof.

Considering any $\phi\in L^{2}(\Gamma_{m,-};|n\cdot v|)$ and $\psi\in L^{2}(\Gamma^{s}_{m,+};|n\cdot v|)$, we have

[TABLE]

where the first equality comes from definition of $\mathcal{P}_{m}$ and the second from the definition of $g$. By multiplying (21) by $f$ and integrating over $\mathscr{D}_{m}^{s}$, we obtain

[TABLE]

where the first equality comes from (21) and second from integration by parts. Comparing the equations above, we have

[TABLE]

It is easy to see that from the definition of $h$ we also have

[TABLE]

By multiplying (22) by $f$ and integrating over $\mathscr{D}_{m}\backslash\mathscr{D}_{m}^{s}$, we obtain

[TABLE]

where the first equality comes from (22) and second from integration by parts. We use notation $\tilde{n}$ as the outer normal direction over $\Gamma^{s}_{m,+}$ and $\Gamma^{s}_{m,-}$ with respect to the domain $\mathscr{D}_{m}\backslash\mathscr{D}_{m}^{s}$, to distinguish it from the outer normal with respect to the domain $\mathscr{D}_{m}^{s}$. In fact, the two instances of “outer normal” have opposite directions when we interpret $\Gamma^{s}_{m,\pm}$ as the boundary of $\mathscr{D}_{m}^{s}$ and the boundary of $\mathscr{D}_{m}\backslash\mathscr{D}_{m}^{s}$. By comparing the equations above, we have

[TABLE]

Noticing that $h=0$ over $\Gamma_{m,+}$ and $\tilde{n}$ is also the outer normal direction over $\Gamma_{m,-}$ when interpreted as the boundary of $\mathscr{D}_{m}$, we have

[TABLE]

where the last equality comes from the definition of $f$ and $\psi$. ∎