Low-Rank Approximation for Multiscale PDEs

2.1. Kinetic equations and asymptotic preserving methods

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')$ is the scattering coefficient that describes the possibility of a photon located at $x$ changing its velocity from $v$ to $v'$, 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 BLP11.

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')$ 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 DG00:

In this case, the photons go through the media that oscillates at a small scale: For example, sunlight passing through a heavy cloud with a large number of small droplets or a 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 DG00.

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 GM01.)

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 it 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 HJL17. 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?

2.2. Elliptic equations and numerical homogenization

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 BL11. 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 All92.

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 EH09 for a discussion of a reduced number of basis functions and E11 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 Hac15; and a Bayesian approach that views the source $f$, and hence the solution $u^\epsilon$, as Gaussian fields Owh15, which further translates to game theory OS19. 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?

3.1. Numerical rank

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 CLLW20a.

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.

3.2. Random sampling in numerical linear algebra

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 numerical PDEs, but it 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 HMT11 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 HMT11.

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.

4.1. Learning basis functions

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)

Offline stage: For $m=1$, $2$, …, $M$, find$$\begin{equation*} \mathsf{G}_m = \left[b_{m,1}\,,b_{m,2}\dots \right]\,, \end{equation*}$$

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{equation} \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} \tag{21}\cssId{texmlid13}{} \end{equation}$$

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

(2)

Online stage: The global solution is$$\begin{equation*} u = \sum _{m=1}^M u_m\chi _m = \sum _{m=1}^M \chi _m\mathsf{G}_{m}{c}_m, \end{equation*}$$

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: