Low-rank approximation for multiscale PDEs

Kinetic equations, which originate from statistical mechanics, describe the evolution of probability density for identical particles in phase space. A model equation, the radiative transfer equation (RTE), characterizes the evolution of photon density. In the steady state, this equation is

where $f(x,v)$ is the light source, and the linear collision operator $\mathsf{S}^{\epsilon}$ describes the interaction of photons with the optical media. The small parameter $\epsilon$ is encoded in this operator.

The operator $\mathsf{S}^{\epsilon}$ defines several distinct regimes. In the optically thick regime, it is defined by

In this case, $k(x,v,v^{\prime})$ is the scattering coefficients that describes the possibility of a photon located at $x$ changing its velocity from $v$ to $v^{\prime}$, and the parameter $\epsilon$ is called the Knudsen number, standing for the ratio of the mean free path to the typical domain length. When the medium is optically thick, the mean free path is small, with $\epsilon\ll 1$. This means the photon particles are scattered fairly often, and the system statistically achieves the equilibrium state, which can itself be characterized mathematically. One example is to observe light in atmosphere, where the average mean free path is about $10$ m, and the observation is conducted at the scale of $10$ km, leading to $\epsilon\sim 10^{-3}$. By performing asymptotic expansion in terms of $\epsilon$, the inhomogeneity in the velocity domain vanishes, and one can show that $u^{\varepsilon}(x,v)$ asymptotically approximates $u^{\ast}(x)$, a function without dependence on $v$ that solves a diffusion equation. We have the following result from [MR2839402].

This result indicates that the homogenized operator as $\epsilon\to 0$ is $\mathcal{L}^{\ast}\propto\nabla_{x}\cdot\left(({1}/{\sigma})\nabla_{x}\right)$. Similar results, when $k(x,v,v^{\prime})$ fails to have the form of $\sigma(x)$ in the anisotropic optical media, are still available, but the explicit form of $\mathcal{L}^{\ast}$ is no longer available.

A second regime of interest for $\mathsf{S}^{\epsilon}$ is one in which the media is highly heterogeneous [MR1760042]:

In this case, the photons go through the media that oscillates at a small scale: For example, sunlight passing through heavy cloud with a large number of small droplets or laser beam passing through crystals. The amplitude of $k$ determines the photon scattering frequency. Since $k$ oscillates rapidly, photons also change rapidly between the high- and low-scattering regimes. On a large scale, the photons can be viewed approximately as scattering with an averaged frequency. A mathematical result is as follows [MR1760042].

In special cases, such as under periodic or random ergodic conditions, the function $\sigma^{\ast}$ can be computed explicitly. (There are also works that investigate the asymptotic limit of the RTE when the system is both highly oscillatory and in diffusion regime; see [MR1878799].)

In both limiting regimes, the limiting equations (8) and (11) can be solved much more efficiently than the original equation (4). The discretization of (4) is constrained strongly by $\epsilon$, due either to stability (as in (5)) or accuracy (as in (9)). By contrast, the solution $u^{\ast}$ varies smoothly, containing no $\epsilon$-scale effects, so can be obtained accurately by applying a discretization with mesh width $O(1)$ to the asymptotic limiting equation. If the latter equation is available, computation of $u^{\ast}$ by this means is the recommended methodology.

Methods for kinetic equations are termed “asymptotic preserving” (AP) if they can relax the requirement that the discretization width $h$ satisfies $h=o(\epsilon)$ yet still capture the asymptotic limits. Many different AP approaches have been proposed. For linear equations, existing AP methods rely on even-odd or micro-macro decomposition. For nonlinear equations, knowledge of the specific forms of the limits is usually required, and this knowledge is built into the solvers [MR3645390]. As mentioned above, these specific forms are often not available, so many AP methods cannot be applied to a large set of multiscale kinetic equations. This observation begs the question: Knowing the existence of the limit, but not its particular form, can we still devise efficient methods for solving kinetic equations?

Another class of multiscale equations that has been investigated deeply is elliptic equations with highly oscillatory coefficient. These problems have the form

where $\epsilon\ll 1$ is the scale on which the media oscillates. (The source term $f$ has no small-scale contribution.)

This equation is a model problem from petroleum engineering where it is crucial to precompute the underground flow before expensive construction of infrastructure takes place [MR2801210]. The problem is typically solved on kilometer-scale domains, but the heterogeneities in the media can scale at centimeters. Certain forms of this equation can be approximated effectively by an equation that can be solved efficiently. Suppose the media coefficient $a^{\epsilon}(x)$ has the form $a(x,x/\epsilon)$, that is, it varies on two scales ($1$ and $\epsilon$), and moreover is periodic with respect to the fast variable (the second argument in $a(x,x/\epsilon)$). Then in the limiting regime as $\epsilon\to 0$, the solution $u^{\epsilon}$ converges to that of a homogenized equation, with the media “smoothed-out,” as described in the following result [MR1185639].

As in the previous section, when a limiting equation can be derived explicitly, the best course for obtaining a useful solution it to solve this equation directly, as the mesh width in the discretization scheme can be much larger than $\epsilon$. See [MR2477579] for a discussion of a reduced number of basis functions and [MR2830582] for computation of the effective media.

However, the validity and the specific form of the effective limit are known only in special cases like the one described in Theorem 3. In other cases, we seek a solver that relies on as little analytical knowledge as possible. An approach known as numerical homogenization has been investigated extensively. This approach is founded on two principles: a discretization scheme independent of $\epsilon$, and a numerical solution scheme that captures the true limiting behavior of the solution on the discrete level. Variants of numerical homogenization include application of the $\mathcal{H}$-matrix, a purely algebraic technique [MR3445676]; and a Bayesian approach that views the source $f$, and hence the solution $u^{\epsilon}$, as Gaussian fields [MR3369060], which further translates to game theory [MR3971243]. All these methods are successful, but they all implicitly rely on properties of the underlying elliptic equation. Can we devise an approach that applies to general problems with oscillatory media that exploits the low-rank property in the solution space, without using analytical structure explicitly?

We consider a bounded linear operator $\mathcal{A}$, which maps $f\in\mathcal{X}$ to a space $\mathcal{Y}$, that is

In the PDE setting, $\mathcal{A}$ is the solution operator that maps the boundary conditions and/or source term $f$ to the solution $u$. The numerical rank of such an operator is defined as follows.

When $\mathcal{A}$ is the PDE solution map, then $\tilde{\mathcal{A}}$ with low rank is also a linear map with a finite dimensionality. It can be viewed as the discrete version (or a matrix of dimension $k_{\tau}(\mathcal{A})$) that approximates $\mathcal{A}$ within $\tau$ accuracy. The definition suggests that if $\tilde{\mathcal{A}}$ can be found, it is optimal in the sense of numerical efficiency. The concept is rather similar to the Kolmogorov $N$-width, defined as follows.

Indeed, the Kolmogorov $N$-width and numerical rank are related by the following result [MR4155236].

For the three examples presented in Section 2, the numerical ranks can be calculated from their limiting equations. For one-dimensional RTE in the diffusion regime, if we denote by $\mathcal{A}^{\epsilon}$ and $\mathcal{A}^{\ast}$ the solution operators of (4) and (8), respectively, then noting that $\mathcal{A}^{\ast}$ can be approximated using $1/\sqrt{\tau}$ grid points to achieve $\tau$ accuracy, when $\epsilon<\tau$, the numerical rank is naturally $k_{\tau}(\mathcal{A}^{\epsilon})\lesssim 1/\sqrt{\tau-\epsilon}$. Without employing the knowledge of the existence of the limit, however, a brute-force discretization naturally requires $O(1/\epsilon\tau^{\alpha+1})$ degrees of freedom: $O(1/\epsilon\tau)$ for the upwind discretization in $x$ and $O(1/\tau^{\alpha})$ for the discretization in $v$, where $\alpha$ depends on the particular numerical integral accuracy. Translating into Green’s-matrix language, this observation means that $\mathsf{G}^{\epsilon}$ is represented by $O(1/\epsilon\tau^{\alpha+1})$ degrees of freedom but its range can be captured by a compressed Green’s matrix $\mathsf{G}^{\ast}$ with just $O(1)$ column vectors.

The same argument applies to the elliptic equation on a two-dimension domain with high oscillations. When second-order linear finite elements are used, with no knowledge of the limiting system, $O(1/\epsilon^{2}\tau)$ degrees of freedom are required, dropping to $O(1/\tau)$ when the the limiting system is known. In other words, the full Green’s matrix $\mathsf{G}^{\epsilon}$ requiring $O(1/\epsilon^{2})$ degrees of freedom can be well-represented using just $O(1)$ column vectors.

In all these cases, the degrees of freedom for a given numerical method are substantially larger than the numerical rank of the problem. Thus, much of the information in these full-blown representations is redundant and compressible. A low rank representation exists and yields a much more economical representation.

Knowing the existence of the low rank structure and finding such a structure are very different goals. The Kolmogorov $N$-width is a concept developed in but has made little impact in numerical PDEs for a simple reason: Traditional PDE solvers require a predetermined set of basis functions, while the Kolmogorov $N$-width looks for “optimal” basis functions. How can an optimal basis be found without first forming the full basis? Translated to linear algebra, this question is about finding the dominant singular vectors in a matrix without forming the whole matrix. Specifically, if $\mathsf{A}\in\mathbb{R}^{m\times n}$ is known to be approximately low rank, meaning that there exists $\mathsf{U}_{r}$, a $m\times r$ matrix with orthonormal columns with $r\ll\min(m,n)$ and

can we find $\mathsf{U}_{r}$ without forming the full matrix $\mathsf{A}$?

In linear algebra, it is well-known that $\mathsf{U}_{r}$ is simply the collection of the first $r$ singular vectors of $\mathsf{A}$. Writing

where $\mathsf{U}=\left[u_{1}\,,u_{2}\,,\dots,u_{n}\right]\in\mathbb{R}^{m\times n}$ and $\mathsf{V}=\left[v_{1}\,,v_{2}\,,\dots,v_{n}\right]\in\mathbb{R}^{n\times n}$ contain the left/right singular vectors and $\Sigma=\mathrm{diag}(\sigma_{1},\sigma_{2},\dotsc,\sigma_{n})$ contains the singular values, then $\mathsf{U}_{r}$ is the first $r$ columns in $\mathsf{U}$.

The standard method for computing the SVD requires $\mathsf{A}$ to be stored and computed with. But the celebrated randomized SVD (rSVD) method [MR2806637] captures the range of a given matrix by means of random sampling of its column space, which requires only computation of matrix-vector products involving $\mathsf{A}$ and random vectors — operations that can be performed without full storage or knowledge of $\mathsf{A}$. Implementation of rSVD is easy and its performance is robust.

The idea behind the algorithm is simple. If matrix $\mathsf{A}\in\mathbb{R}^{m\times n}$ has approximate low rank $r\ll\min\{m,n\}$, the matrix maps an $n$-dimensional sphere to an $m$-dimensional ellipsoid that is “thin:” $r$ of its axes are significantly larger than the rest. With high probability, vectors that are randomly sampled on the $n$-dimensional sphere are mapped to vectors that lie mostly in a $r$-dimensional subspace of $\mathbb{R}^{m}$ — the range of $\mathsf{A}$. An approximation to $\mathsf{A}_{r}$ can be obtained by projecting onto this subspace.

A precise statement of the performance of randomized SVD is as follows [MR2806637].

The result reconstructs the range of $\mathsf{A}$ in a nearly optimal way. It is optimal in efficiency because to capture a rank-$r$ matrix, only $r+p$ matrix-vector products involving $\mathsf{A}$ are required for the calculation of $\mathsf{Y}$, and the oversampling parameter $p$ is typically quite modest. ($p=5$ is a typical value.) The result is nearly optimal in accuracy as well. The error bound relies only on $\sigma_{r+1}$, which is expected to be smaller than $\sigma_{1}$. The decay profile of singular values do not affect the approximation accuracy.

If a low-rank approximation to the matrix $\mathsf{A}$ is required, and not just an approximation of its range, another step involving multiplications with its transpose is needed. The full method is shown in Algorithm 1.

In standard domain decomposition, the local discretized Green’s matrix $\mathsf{G}_{m}$ is assembled from a “full” collection of basis functions in the patch $\mathcal{K}_{m}$. The global solution to (18), confined to each $\mathcal{K}_{m}$, is a linear combination of the columns $\mathsf{G}_{m}$. The coefficients of these combinations are chosen so that that the continuity conditions across patches and the exterior boundary condition are all satisfied. The complete process can be outlined as follows.

1. (1)

   Offline stage: For $m=1,2,\dots,M$, find

   $$\mathsf{G}_{m}=\left[b_{m,1}\,,b_{m,2}\dots\right]\,,$$

   where each local function $b_{m,n}$ is a solution to (18) restricted to the subdomain $\mathcal{K}_{m}$, with fine grid $h\ll\epsilon$ and delta-function boundary conditions. That is,

   $$\begin{cases}\mathcal{L}^{\epsilon}b_{m,n}=0\,,\quad&x\in\mathcal{K}_{m}\\ b_{m,n}=\delta_{m,n}\,,\quad&x\in\partial\mathcal{K}_{m}\,,\end{cases}$$

   (21)

   where $\delta_{m,n}$ is the Kronecker delta that singles out the $n$-th grid point on the boundary $\partial\mathcal{K}_{m}$.

2. (2)

   Online stage: The global solution is

   $$u=\sum_{m=1}^{M}u_{m}\chi_{m}=\sum_{m=1}^{M}\chi_{m}\mathsf{G}_{m}{c}_{m},$$

   with the support of each $u_{m}=\mathsf{G}_{m}{c}_{m}$ confined to $\mathcal{K}_{m}$, where $c_{m}$ is a vector of coefficients determined by the boundary conditions $\phi$ and continuity conditions across the patches.

The complete basis represented by $\mathsf{G}_{m}$ has a low-rank structure that can be revealed using randomized SVD. Instead of using delta functions as the boundary conditions, we propose to obtain basis functions by setting random values on $\partial\mathcal{K}_{m}$, as follows:

where $\omega_{m,n}$ is defined to have a random value drawn i.i.d. from a normal distribution at each grid point in $\partial\mathcal{K}_{m}$. Denoting $\mathsf{G}_{m}^{\mathrm{r}}=\{r_{m,1}\,,r_{m,2}\,,\cdots\}$, we have from linearity of the equation that

where $\Omega$ is a random i.i.d. matrix with entries $\omega_{m,n}$. This $\mathsf{G}_{m}^{\mathrm{r}}$ is used in the online stage, as an accurate surrogate of $\mathsf{G}_{m}$, see Algorithm 2.

Although we do not apply full-blown rSVD here, the homogenizable and low-rank property of the local solution space implies that $\mathsf{G}_{m}^{\mathrm{r}}$ and $\mathsf{G}_{m}$ share similar range with the number of basis functions $k_{m}$ in $\mathsf{G}_{m}^{\mathrm{r}}$ being much smaller than $n_{m}$, the number of grid points on $\partial\mathcal{K}_{m}$, and independent of $\epsilon$. In Figure 3 we plot the angles between $\mathsf{G}_{m}^{\mathrm{r}}$ and $\mathsf{G}_{m}$ for two of the equations discussed in Section 2. In both cases, and for small $\epsilon$, the approximated Green’s matrix quickly recovers the true Green’s matrix as the number of samples $k_{m}$ increases, and thus captures the local solution space. We should note that if $\mathsf{G}_{m}$ does not have low-rank structure, in the sense that $k_{m}\sim n_{m}$, then solving (22) would be equally expensive as solving (21), hence the random sampling technique does not gain any computational efficiency when the system is not homogenizable.

In Figure 4 we showcase the basis functions on a patch for elliptic equation with media coefficient $a(x_{1},x_{2})=1+1000\,\mathbf{1}_{S}(x_{1},x_{2})$, with $S=\{(x_{1},x_{2})\in[0,1]^{2}:(x_{1}\cos(100\sqrt{(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}}))\leq x_{2}-0.5\}$. For this non-conventional media without any periodic structure, the traditional multiscale methods are no longer valid, but our method still quickly captures the optimal basis.

For particular boundary conditions $\phi$, the global solution is assembled from the local basis functions in the online stage. Two numerical examples are shown in Figure 5. In both examples, there is little visible difference between the reference solution and the approximated one computed from the reduced random basis. Only 8.3% and 62.5%, respectively, of the degrees of freedom required by the full basis are needed to represent these solutions using a random basis; see details in [MR4155236].

Since we do not have access to the full set of basis functions, the condition that $u$ defined in (23) is continuous across subdomains can be satisfied only in a least-squares sense; see [MR4155236] for details.

Our second approach for exploiting the low-rank property in multiscale computations is based on Schwarz iteration. The Schwarz method is a standard iteration algorithm within the domain decomposition framework, in which boundary-value problems are solved on the patches, with neighboring patches subsequently exchanging information and re-solving until consistency is attained. The exchange of boundary information between neighboring patches is known as the boundary-to-boundary (BtB) map. The map has an exploitable low-rank property.

To develop the approach, we write the solution of (18) as

where the partition-of-unity functions $\chi_{m}$ are defined in (20). The solution $u^{\epsilon}_{m}$ on patch $m$ is uniquely determined by $f_{m}$, its local boundary condition, according to the equation

The Schwarz method starts by initial guesses to the local boundary conditions $f_{m}=f_{m}^{0}$, then on iteration $t$, it solves the subproblems (25) with $f_{m}=f_{m}^{t}$ to obtain all local solutions $u_{m}^{\epsilon}$. By confining $u_{m}^{\epsilon}$ on the boundaries of adjacent patches, one updates the boundary conditions for surrounding patches:

Here $\mathcal{S}_{m}$ denote the solution to (25) confined in the interior of $\mathcal{K}_{m}$, and $\mathcal{P}_{m}$ takes the trace of the solution on the neighboring boundaries $\mathcal{K}_{m}\cap\partial\mathcal{K}_{n}$ for $n\in\mathcal{I}_{m}$, for the updated boundary condition.

Define the BtB map by $\mathcal{A}_{m}:=\mathcal{P}_{m}\circ\mathcal{S}_{m}$, and define $\mathcal{A}$ and $f^{t}$ to be the aggregation of $\mathcal{A}_{m}$ and $f_{m}^{t}$, respectively, over $i=1,2,\dotsc,M$. We can then write the updating procedure as

The overall method is summarized in Algorithm 3.

Most of the computation in the Schwarz method during the iteration comes from solving the boundary-value PDEs on the patches, to implement the map $\mathcal{A}$. Since the PDE is homogenizable, the solution space on each patch is approximately low rank, and the map $\mathcal{A}$ can be expected to inherit this property. If we can “learn” this operator in an offline stage, and simply apply a low-rank approximation repeatedly in the online stage, the online part of Algorithm 3 can be made much more efficient. In our approach, Algorithm 1 is used to compress the map $\mathcal{A}$.

This approach is quite different from the one described in Section 4.1, in the sense that it is not only the range of the solution space, but the whole operator that is being approximated. To apply Algorithm 1, we need to define the “adjoint operator” for $\mathcal{A}$ on the PDE level. This operator is composed of the adjoints $\mathcal{S}_{m}^{\ast}$ for the local solution operators $\mathcal{S}_{m}$ of (25) on each domain $\mathcal{K}_{m}$. The form of $\mathcal{S}_{m}^{\ast}$ is specific to the PDE; we use the elliptic equation as an example. Defining $\mathcal{L}^{\epsilon}=\nabla\cdot\left(a(x,x/\epsilon)\nabla\right)$, $\mathcal{S}_{m}^{\ast}$ is defined in the following result.

We also describe calculation of the adjoint operator $\mathcal{S}_{m}^{\ast}$ for the RTE (4).

The specific form of the adjoint operator $\mathcal{S}_{m}^{\ast}$ allows us to adapt the randomized SVD algorithm to compress the confined solution map $\mathcal{S}_{m}$; see Algorithm 4. This method requires only $k$ solves of local PDE (25) and sourced adjoint PDE (28) (or (30)), together with a QR factorization and SVD of relatively small matrices. The overall low-rank Schwarz iteration is then summarized in Algorithm 5.

In Figure 6 we present the confined solution operator $\mathcal{S}_{m}$ and its low-rank approximation $\mathcal{S}_{m}^{\mathrm{r}}$. For the radiative transfer equation, the map is a 2880-by-40 matrix, with the size of each patch being 0.2$\times$[-1,1], with $\Delta x=1/360$ and $\Delta v=1/40$. The random sampling procedure reconstructs it with just 6 samples. For the elliptic equation, the map is a 1600-by-160 matrix, and the size of each patch is 1$\times$[0,1] with $\Delta x=1/40$. The random sampling approximates it well with 60 samples. The compression rates for these examples are thus 6.7 and 2.7, respectively. See details in [MR4252068].

In Figure 7, we show numerical examples for the global solutions of two problems obtained from the approach of this section. The reference solution (obtained with a fine mesh) is well captured by the approximation that uses the low-rank BtB map as the surrogate in the Schwarz iteration. These two cases use just $15\%$ and $43\%$, respectively, of the number of local solves needed to capture the BtB map at fine scale. While the relative error of the reduced Schwarz method decays as fast as the standard Schwarz iteration, as shown in Figure 8, the cost is much reduced. See Table 1 for a comparison of computation times.