Model input and output dimension reduction using Karhunen Lo\`{e}ve expansions with application to biotransport

TL;DR

This paper investigates dimension reduction techniques for biotransport modeling in tumors, demonstrating that pressure fields can be accurately represented with few KL expansion terms despite high-dimensional input parameters, enabling faster simulations.

Contribution

It introduces a reduced-order modeling approach using KL expansions to efficiently simulate tumor biotransport with uncertain heterogeneous properties.

Findings

Pressure fields are less sensitive to high-order KL terms of input.

Few KL terms suffice to accurately represent pressure fields.

High-dimensional input representations are needed for permeability, but not for pressure.

Abstract

We consider biotransport in tumors with uncertain heterogeneous material properties. Specifically, we focus on the elliptic partial differential equation (PDE) modeling the pressure field inside the tumor. The permeability field is modeled as a log-Gaussian random field with a prespecified covariance function. We numerically explore dimension reduction of the input parameter and model output. Truncated Karhunen--Lo\`{e}ve (KL) expansions are used to decompose the log-permeability field, as well as the resulting random pressure field. We find that although very high-dimensional representations are needed to accurately represent the permeability field, especially in presence of small correlation lengths, the pressure field is not very sensitive to high-order KL terms of the input parameter. Moreover, we find that the pressure field itself can be represented accurately using a KL expansion with a small number of terms. These observations are used to guide a reduced-order modeling approach to accelerate computational studies of biotransport in tumors.

Figures

Equations

$$\underset{\mathbf{h}\to 0}{\lim }E[{\left(Z(\mathbf{x}+\mathbf{h},\cdot )-Z(\mathbf{x},\cdot )\right)}^2]=0,\text{for all }\mathbf{x}\in X.$$

$$c(\mathbf{x},\mathbf{y})$$

$$\rho (\mathbf{x},\mathbf{y})$$

$$[Cu](\mathbf{x})={\int }_{X}c(\mathbf{x},\mathbf{y})u(\mathbf{y})d\mathbf{y},u\in L^2(X).$$

$${\int }_{X}c(\cdot ,\mathbf{y})e_i(\mathbf{y})d\mathbf{y}={\lambda }_i{e}_i(\cdot ),i=1,2,\dots .$$

$$Z(\mathbf{x},\omega )=\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\xi }_i(\omega )e_i(\mathbf{x}),$$

$${\xi }_i(\omega )=\frac{1}{\sqrt{{\lambda }_i}}{\int }_{X}Z(\mathbf{x},\omega )e_i(\mathbf{x})d\mathbf{x}.$$

$$U(\mathbf{x},\mathbf{\xi })$$

$$\bar{U}(\mathbf{x})$$

$$u_i(\mathbf{\xi })=\frac{1}{\sqrt{{\lambda }_i}}{\int }_X(U(\mathbf{x},\mathbf{\xi })-\bar{U}(\mathbf{x}))e_i(\mathbf{x})d\mathbf{x},$$

$${\bar{U}}_m=\frac{1}{N}\sum _{j=1}^{N}U({\mathbf{x}}_m,{\mathbf{\xi }}^j),m\in \{1,\dots ,N_{\text{quad}}\}.$$

$$u_c({\mathbf{x}}_m,{\mathbf{\xi }}^k)=U({\mathbf{x}}_m,{\mathbf{\xi }}^k)-{\bar{U}}_m,k\in \{1,\dots ,N\},m\in \{1,\dots ,N_{\text{quad}}\}.$$

$$K_{lm}=\frac{1}{N-1}\sum _{k=1}^N{u}_c({\mathbf{x}}_l,{\mathbf{\xi }}^k)u_c({\mathbf{x}}_m,{\mathbf{\xi }}^k),l,m\in \{1,\dots ,N_{\text{quad}}\}.$$

$${\mathbf{W}}^{1/2}\mathbf{K}{\mathbf{W}}^{1/2}{\mathbf{v}}_i={\lambda }_i{\mathbf{v}}_i,i\in \{1,\dots ,N_{\text{quad}}\}.$$

$$u_i({\mathbf{\xi }}^k)=\frac{1}{\sqrt{{\lambda }_i}}\sum _{m=1}^{N_{\text{quad}}}{w}_m{u}_c({\mathbf{x}}_m,{\mathbf{\xi }}^k)e_i^m,i\in \{1,\dots ,N_{kl}\},k\in \{1,\dots ,N\}.$$

$$-\frac{d}{dx}\left(\kappa (x,\omega )\frac{dp(x,\omega )}{dx}\right)$$

$$p(-1,\omega )$$

$$p(1,\omega )$$

$$c_{Z}(x,y)=\exp\Big{\{}\!\!-\!\frac{|x-y|}{\ell}\Big{\}}.$$

$$a(x,\omega )=a_0(x)+\sigma Z(x,\omega ),$$

$$a_{N_{kl}^a}(x,\omega )=a_0+\sigma \sum _{i=1}^{N_{kl}^a}\sqrt{{\lambda }_i}{\xi }_i(\omega )e_i(x),$$

$$\mathbf{\xi }={\left[\begin{array}{cccc}{\xi }_1& {\xi }_2& \cdots & {\xi }_{N_{kl}^a}\end{array}\right]}^T$$

$$p(x,\mathbf{\xi })=\bar{p}(x)+\sum _{j=1}^{\infty }\sqrt{{\lambda }_j(C_p)}{p}_j(\mathbf{\xi })v_j(x)$$

$$p_j(\mathbf{\xi })=\frac{1}{\sqrt{{\lambda }_j(C_p)}}{\int }_D(p(x,\mathbf{\xi })-\bar{p}(x))v_j(x)dx,j=1,2,\dots ,$$

$$\mathrm{Var}\left[\sum _{j=1}^{N_{kl}^p}\sqrt{{\lambda }_j(C_p)}{p}_j(\mathbf{\xi })v_j(x)\right]=\sum _{j=1}^{N_{kl}^p}{\lambda }_j(C_p)v_j(x{)}^2.$$

$$\frac{\partial }{\partial r}\left(\frac{\kappa r}{\mu }\frac{\partial p}{\partial r}\right)+\frac{1}{r}\frac{\partial }{\partial \theta }\left(\frac{\kappa }{\mu }\frac{\partial p}{\partial \theta }\right)=0.$$

$$\{\begin{array}{rlr}p& =0,& r=R_{\text{tumor}},\\ \frac{\partial p}{\partial r}& =\frac{-Q\mu }{2\pi R_{\text{needle}}\kappa },& r=R_{\text{needle}}.\end{array}$$

$$a(\mathbf{x},\omega )=a_0(\mathbf{x})+{\sigma }_{a}Z(\mathbf{x},\omega ),\mathbf{x}\in D,\omega \in \Omega .$$

$$a(\mathbf{x},\omega )\approx a_0(\mathbf{x})+{\sigma }_a\sum _{i=1}^{N_{kl}^a}{\sqrt{\lambda }}_i{\xi }_i(\omega )e_i(\mathbf{x}),$$

$$p(\mathbf{x},\mathbf{\xi })=\bar{p}(\mathbf{x})+\sum _{j=1}^{N_{kl}^p}\sqrt{{\lambda }_j(C_p)}{p}_j(\mathbf{\xi })v_j(\mathbf{x}),$$

Full text

1 INTRODUCTION

We focus on modeling biotransport processes in tumors with uncertain heterogeneous material properties. An improved understanding of these processes can provide vital insight for agent delivery in cancer treatment [1, 2]. Biotransport processes in tumors can be modeled as flows in heterogeneous porous media. Equations governing biotransport consist of an elliptic PDE describing the pressure distribution and a hyperbolic PDE that describes agent (e.g., drug) delivery in porous media [3]. The uncertain tumor material properties can be modeled as random fields, which are then incorporated as coefficient functions in the governing PDEs.

In the present work, the uncertain permeability field is modeled as a log-Gaussian random field. Our aim is to efficiently simulate the uncertain pressure field. We use KL expansions [4, 5] to represent the random log-permeability field. The use of KL expansions for representing random field parameters in mathematical models has been a common modeling approach in the uncertainty quantification community [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The typical approach is to use a truncated KL expansion with enough term to ensure the average variance of the parameter field is sufficiently captured. That is, the truncation of the KL expansion is performed a priori and without taking the response of the governing equations to the random field coefficients in mind. We take a goal-oriented point of view: instead of relying on a truncated KL expansion of the log-permeability field that is computed independently of the governing PDE, we seek to retain only the KL terms that the PDE solution operator is sensitive to. This goal-oriented strategy can lead to significant input parameter dimension reduction, especially for input fields with small correlation lengths. The PDE solution—the pressure field—itself can also be represented via a truncated KL expansion. We observe that a low-rank representation of the pressure field is often afforded by a truncated KL expansion with a small number of terms. The latter is a consequence of the (often) rapid decay of the eignvalues of the model output covariance operator. Our approach guides an input and output dimension reduction strategy: a low-rank representation of the pressure field can be computed in a low-dimensional parameter space. We mention that a preliminary version of this work was presented in the conference paper [16].

This article is structured as follows. In Section 2, we recall the requisite background material on random fields and their KL expansion. In that section, we also outline a computational strategy for computing KL expansions for random fields with or without a prespefied covariance function. In Section 3, we use a model elliptic PDE in one space dimension to illustrate the components of the proposed approach. Then, in Section 4, we focus on a biotransport application problem. We present numerical results illustrating the merits of the proposed strategy. Concluding remarks are provided in Section 5.

2 BACKGROUND ON RANDOM FIELDS

Let $(\Omega,\mathcal{F},P)$ be a probability space, where $\Omega$ is a sample space, $\mathcal{F}$ is an appropriate $\sigma$-algebra, and $P$ is a probability measure. Let $X\subset\mathbb{R}^{d}$, with $d=1,2$, or $3$, be a compact set. Let $Z:X\times\Omega\rightarrow\mathbb{R}$ be a stochastic process [17]. From a modeling standpoint, $Z(\mathbf{{x}},\omega)$ can be used to represent uncertain parameters fields in mathematical models.

A stochastic process is called centered if $E[Z(\mathbf{{x}},\cdot)]=0$ for all $\mathbf{{x}}\in X$, where $E[Z(\mathbf{{x}},\cdot)]=\int_{\Omega}Z(\mathbf{{x}},\omega)\,P(d\omega)$. A process $Z$ is called mean square continuous if

[TABLE]

The covariance function $c:X\times X\rightarrow\mathbb{R}$ and the corresponding correlation function of a stochastic process $Z$ are, respectively, given by

[TABLE]

We also recall the definition of the covariance operator of a stochastic process $Z(\mathbf{{x}},\omega)$, which is given by

[TABLE]

Karhunen–Loève expansion. Let $Z:X\times\Omega\rightarrow\mathbb{R}$ be a centered mean-square continuous stochastic process, and let $\{e_{i}\}_{i=1}^{\infty}$ be the orthonormal basis of eigenvectors of its covariance operator with corresponding (non-negative) eigenvalues $\{\lambda_{i}\}_{i=1}^{\infty}$:

[TABLE]

The process $Z(\mathbf{{x}},\omega)$ can be represented via its KL expansion [4, 5, 18, 19]:

[TABLE]

where $\xi_{i}$ are centered mutually uncorrelated random variables with unit variance and are defined by

[TABLE]

The convergence of the series (3) is uniform in $X$, and is mean square in $\Omega$ [4]. Moreover, if $Z(\mathbf{{x}},\omega)$ is a Gaussian process, convergence of the series  (3) is almost sure for each $\mathbf{{x}}\in X$; see [4, p. 485] for further details.

Numerical computation of KL expansion. To compute the KL expansion of a stochastic process the eigenvalue problem (2) must be solved first. In the present work, we follow Nyström’s approach [20], which involves discretizing the generalized eigenvalue problem using quadrature. We describe the steps for computing KL expansions below. Further details on numerical methods for computing KL expansions can be found in [21].

When modeling random field coefficients in models, one often has access to a prespecified covariance function. On the other hand, when computing KL expansion of a random field output of a mathematical model we only have access to realizations of the model output. Let $U$ denote the random field output of model governed by PDEs. In practice, often the model uncertainties are parameterized using a random vector $\bm{\xi}$, in which case the random field output $U=U(\mathbf{{x}},\bm{\xi})$ can be computed for specific realizations of $\bm{\xi}$. To compute the truncated KL expansion,

[TABLE]

the covariance function of $U$ needs to be approximated via sampling, resulting in an approximate covariance operator $C$ for the process. Then, the generalized eigenvalue problem will be solved using this approximate covariance operator to find (approximations to) $\lambda_{i}$ and $e_{i}$, $i=1,\ldots,{N_{kl}}$. In practice, the dominant KL terms can be captured reliably, with a modest sample size, as discussed in our numerical results. We summarize the steps required for computing truncated KL expansion of the random process $U(\mathbf{{x}},\bm{\xi})$ in Algorithm 1. In what follows, we refer to the coefficients $u_{i}$ in (4) as the KL modes.

3 MODEL 1D ELLIPTIC EQUATION WITH RANDOM COEFFICIENT FUNCTION

We let $(\Omega,\mathcal{F},P)$ be a probability space, and for $\omega\in\Omega$, consider the following model elliptic boundary value problem:

[TABLE]

In the following numerical experiments, the right hand side function is given by $f(x)=\cos(\pi x)+\sin(2\pi x)$. We model the coefficient function $\kappa(x,\omega)$ as a log-Gaussian random field as follows. Let $Z(x,\omega)$ be a centered Gaussian process with covariance function,

[TABLE]

We set the correlation length $\ell$ of the process to $\ell=1/4$. Then, we set $\kappa(x,\omega)=\exp(a(x,\omega))$ with

[TABLE]

where $a_{0}$ and $\sigma^{2}$ are the pointwise mean and variance of $a(x,\omega)$, respectively; we choose these parameters such that the pointwise mean and standard deviation of $\kappa(x,\omega)$ are $m=0.1$ and $s=.07$, respectively. Accordingly, we let $\sigma^{2}=\log(1+s^{2}/m^{2})$ and a_{0}\equiv\log\big{(}m/\sqrt{1+s^{2}/m^{2}}\big{)}. 111 We have used the well-known formulas relating the mean and variance of a log-normal random variable $Y=\exp(a_{0}+\sigma X)$, where $X$ is standard normal, to $a_{0}$ and $\sigma^{2}$. We consider the weak formulation of the problem (5), and use the continuous Galerkin finite element method, with linear basis functions, to solve the problem numerically.

We use a truncated KL expansion for $a(x,\omega)$,

[TABLE]

where $(\lambda_{i},e_{i})$ are eigenpairs of the covariance operator of $Z(x,\omega)$. Due to the Gaussianity of the process, $\xi_{i}$ are independent standard normal random variables. Note that the random vector

[TABLE]

completely parameterizes the uncertainty in the problem (5), and its solution $p({x},\omega)=p(x,\bm{\xi}(\omega))$.

As a first illustration, we consider a fixed realization of $a(x,\omega)$ as ${N_{kl}^{a}}$ increases in Figure 1 (left). Note that sufficiently large ${N_{kl}^{a}}$ is needed to capture the fluctuations of the random field. On the other hand, the corresponding PDE solution is less sensitive to the higher-order KL terms of the parameter, as seen in Figure 1 (right). This behavior is consistent with the analysis in [22], where a global sensitivity analysis formalism is used to quantify the impact of the KL terms of the log-coefficient, in an elliptic PDE, on variability in solution of the PDE.

Next, we study the properties of the PDE solution $p(x,\bm{\xi})$. We depict the correlation function of $p(x,\bm{\xi})$, approximated via Monte Carlo sampling (with $10^{4}$ samples), in Figure 2 (top left). This indicates strong correlations in the output field $p(x,\bm{\xi})$. In Figure 2 (top right) we compare the (normalized) eigenvalues of the covariance operators for $p(x,\bm{\xi})$ and $a(x,\omega)$; we note a much faster spectral decay for the output covariance operator. The latter indicates that a KL expansion with a small number of terms can be used to approximate $p(x,\bm{\xi})$ reasonably well. We study this by considering the KL expansion

[TABLE]

of $p(x,\bm{\xi})$, where $(\lambda_{j}(C_{p}),v_{j})$ are the eigenpairs of covariance operator $C_{p}$ of $p$, computed numerically using Algorithm 1, $p_{j}$ are given by

[TABLE]

and $\bar{p}(x)$ is the mean of $p(x,\bm{\xi})$.

To quantify the impact of truncating the KL expansion of $p(x,\bm{\xi})$ on its approximation properties, we study the pointwise variance $\mathrm{Var}[p(x,\bm{\xi})]$ with different choices of ${N_{kl}^{p}}$. Note also that

[TABLE]

The results in Figure 2 (bottom left), indicate that pointwise variance of $p(x,\bm{\xi})$ can be approximated well with a small ${N_{kl}^{p}}$. Note that the finite-element grid resolution used to solve the PDE also affects the accuracy the KL expansion of the output. In Figure 2 (bottom right) we perform a grid refinement study as we compute the pointwise variance of the process, where we fix ${N_{kl}^{p}}={N_{kl}^{a}}=10$. For the present problem using about 50 grid points seems to be sufficient to resolve the pointwise variance. More broadly, one needs a sufficiently fine computational grid to ensure the dominant eigenpairs of the covariance operator are resolved with sufficient accuracy. The grid resolution issues become more consequential in problems in two or three space dimensions, as the dimension of the discretized eigenvalue problem can become quite large.

We next study input parameter and output dimension reduction in Figure 3 where we show typical realizations of $p(x,\bm{\xi})$, for a small ${N_{kl}^{a}}$ (top row) and a relatively large ${N_{kl}^{a}}$ (bottom row) for various choices of ${N_{kl}^{p}}$ (output dimension).

The numerical experiments in this section lead to the following observations: (i) it is possible to reduce parameter dimension by focusing on KL terms of the parameter that the PDE solution operator is most sensitive to; and (ii) it is possible to reduce output dimension by focusing on the dominant KL terms of the output. In the next section, we explore these notions systematically, in a more challenging problem, involving biotransport in tumors.

4 APPLICATION TO BIOTRANSPORT IN TUMORS

Governing equations and numerical setup. In this section, we study the pressure field in a tumor when a single needle injection occurs at the tumor center. A 2D model in a polar coordinate system is used to analyze the flow field. Consider the mass conservation law and Darcy’s law for steady incompressible flows in a 2D domain, $D=\{(r,\theta):R_{\text{needle}}<r<R_{\text{tumor}},\,0<\theta<2\pi\}$,

[TABLE]

Here $p$ is the pressure, $\kappa$ is the permeability, $\mu$ is the fluid dynamic viscosity, $r$ is the radial distance from a fixed origin, $\theta$ is the polar angle, $R_{\text{tumor}}$ is the radius of the tumor, and $R_{\text{needle}}$ is the radius of the needle used to inject nanofluid into the tumor.

The boundary conditions for the pressure equation are specified as follows:

[TABLE]

Herein, $Q$ is the volume flow rate per unit length. Periodic boundary conditions are enforced in the $\theta$ direction. In this study, $R_{\text{needle}}$ and $R_{\text{tumor}}$ are set to $0.25~{}mm$ and $5~{}mm$, respectively, $Q$ is $1~{}mm^{2}/min$, and $\mu$ is $8.9\times 10^{-4}~{}Pa\cdot s$.

Uncertainties in permeability field. As before, let $(\Omega,\mathcal{F},P)$ be a probability space. Following [23], the permeability field $\kappa$ is modeled by a log-Gaussian random field, and its mode is set to $0.5~{}md$, where $md$ stands for millidarcy. We assume that the log-permeability, a(\mathbf{{x}},\omega)=\log\big{(}\kappa(\mathbf{{x}},\omega)), is given by

[TABLE]

Here $a_{0}$ is the pointwise mean of the process, $\sigma^{2}_{a}$ is the pointwise variance, and $Z$ is a centered Gaussian process with unit pointwise variance for every $\mathbf{{x}}\in D$. In this study, $\sigma_{a}^{2}$ is set to $0.25$, and $a_{0}$ is calculated from the definition of the mode of $\kappa$ as $a_{0}=\ln(0.5)+\sigma_{a}^{2}$. The covariance function of $Z$ is expressed as $c_{Z}(\mathbf{{x}},\mathbf{{y}})=\exp\left\{-\frac{1}{\ell}\|\mathbf{{x}}-\mathbf{{y}}\|_{1}\right\}$, $\mathbf{{x}},\mathbf{{y}}\in D$, where $\ell>0$ is the correlation length. As before, the (Gaussian) log-permeability field can be expressed with a truncated KL expansion,

[TABLE]

where $\lambda_{i}$ and $e_{i}$ are eigenpairs of the covariance operator of $Z$, and $\xi_{i}$ are independent standard normal random variables.

Uncertainty properties of the pressure field. Whereas the governing equation is more complex than the previously considered 1D problem, it is still an elliptic PDE and hence we observe similar behavior in terms of potential for dimension reduction. In the present example, we focus on pressure field over regions around the injection site. Specifically, we consider annular regions with inner boundary given by the inner boundary of the domain and the outer boundary specified by circles of radius $R_{\text{qoi}}=1$ mm, $R_{\text{qoi}}=2$ mm, or $R_{\text{qoi}}=3$ mm. Three sets of realizations of the permeability field (in the entire domain) and the corresponding model output (in the annular domain with $R_{\text{qoi}}=3$ mm) are presented in Figure 4. We observe that although the permeability field realizations exhibit complicated features, the fluctuations in the pressure field are mild.

We denote the covariance operator of the log-permeability field by $C_{a}$ and that of the pressure field by $C_{p}$. In Figure 5 (top), we report the (normalized) eigenvalues of $C_{a}$ and those of $C_{p}$, corresponding to $R_{\text{qoi}}=1$ mm, $R_{\text{qoi}}=2$ mm, and $R_{\text{qoi}}=3$ mm. First, we note that the eigenvalues of $C_{p}$ exhibit a far more rapid decay as compared to that of $C_{a}$. Moreover, as the size of region of interest decreases, the spectral decay of $C_{p}$ becomes faster. In Figure 5 (middle), we examine the spectral decay of $C_{p}$, as the correlation length of the log-permeability increases; for this test we used $R_{\text{qoi}}=3$ mm. As expected, increasing the complexity of the input parameter, by using smaller correlation lengths, leads to slower spectral decay for $C_{p}$. However, we find that even with $\ell=0.5$ mm, about $96\%$ of average output variance is captured by the first 20 KL modes of the output. Finally, Figure 5 (bottom) summarizes the effect of input and output dimensions on capturing the average variance of the output (with $R_{\text{qoi}}=3$ mm, and input parameter correlation length $\ell=1$ mm). Note that the variance of $p(\mathbf{{x}},\bm{\xi})$, restricted to the region of interest, is computed by $\sum_{k=1}^{N_{kl}^{p}}\lambda_{k}(C_{p})$. We note that the average variance can be approximated with reasonable accuracy with small ${N_{kl}^{a}}$ and ${N_{kl}^{p}}$.

We also examine the average relative $L^{2}$ error of the truncated KL representation of the output (with $R_{\text{qoi}}=3$ mm) as ${N_{kl}^{a}}$ and ${N_{kl}^{p}}$ increase, for input (i.e., permeability) fields with different correlation lengths; the results are reported in top and bottom panels of Figure 6, respectively. For the figure in top, we used the KL expansion of input with $2{,}000$ terms as a reference true log-permeability field. For the figure at the bottom, we compute the relative $L^{2}$ error of the output KL representation with the PDE solution restricted to the region of interest. We note that when the input dimension is fixed, the average relative error of the output KL expansion decreases very fast when the number of output KL modes increase, and is not very sensitive to input parameter correlation length. On the other hand, for small correlation lengths, there is a notable increase in the number of input KL modes needed to represent the output accurately.

Insights into reduced-order modeling. From previous analysis, we observed that the spectrum of the output covariance operator $C_{p}$ decays very fast, even when the correlation length $\ell$ is small. This indicates that the output, i.e., the pressure field, can be effectively approximated by a truncated KL expansion as

[TABLE]

with a small number ${N_{kl}^{p}}$ of KL terms. The importance of this approximation is that it decouples the spatial (i.e., $\mathbf{{x}}$) dimensions and those of the random variable $\bm{\xi}(\omega)$. If a surrogate model, such as a polynomial chaos expansion (PCE) [5, 18] is directly used to approximate the pressure field for the purpose of uncertainty quantification [23], the PCE needs to be built for each spatial point on the computational mesh. Instead, if the approximation in (15) is used, PCE (or any other suitable surrogate model) only needs to be constructed for each KL mode $p_{j}$, $j=1,\ldots,{N_{kl}^{p}}$.

We evaluate the performance of the reduced-order model (ROM), i.e., (15), on recovering the probability density functions (PDFs) of pressures at different locations in the flow field. The KL expansion of output is computed using Algorithm 1. As shown in [16], a modest sample size $N$ can be used to capture the dominant modes of output KL expansion. Here, to ensure accuracy, we use $N=1{,}000$ samples. The input dimension is fixed at ${N_{kl}^{a}}=150$ in following tests. In Figure 7, we present four points, namely, $P_{1}$, $P_{2}$, $P_{3}$ and $P_{4}$, on the mesh where PDFs are constructed. The contour stands for relative standard deviation (RSD) of the pressure field. Two correlation lengths, namely, $\ell=0.5$ mm and $\ell=3$ mm, are tested. In both cases, the region of the quantity of interest has a outer radius $R_{\text{qoi}}$ of $3$ mm.

In Figure 8, we present the PDFs from the full model, which is the numerical solution of the governing equation (12), and those from the ROMs with the first $n$ KL terms, where $n\in\{1,5,10,20,40\}$, when $\ell$ is $0.5$ mm. We observe that although the permeability field in this case is very complex, almost all ROMs with the first 10 KL terms can reasonably recover the PDFs constructed from the corresponding full models. With the first 40 KL terms, ROMs can recover the PDFs constructed from full models with negligible discrepancy on all the four points studied. When the correlation length becomes larger, e.g., $\ell=3$ mm, the effectiveness of ROMs becomes more apparent than that with small correlation lengths. From Figure 9, we can clearly see that when $\ell$ equals to $3$ mm, at the point $P_{1}$, the ROM with only the first KL mode can almost recover the PDF constructed from the full model; at the point $P_{4}$, the ROM with the first five KL modes can capture almost all the features in the PDF.

As in the 1D model problem examined earlier, we can achieve a substantial output dimension reduction by using the ROM given by it KL expansion. However, the case for the input dimension reduction is less clear in this case. When the correlation length $\ell$ is small, a very high input dimension is needed to capture most of the variance in the permeability field. However, as before, the PDE solution is still not very sensitive to high-order KL modes; for instance, even with $\ell=0.5$ mm, we note that the average relative error falls below $5\%$ with an ${N_{kl}^{a}}$ of around $50$ (see Figure 6). However, if further accuracy is required, more input KL terms need to be retained. A question arises: is there a way to find a subset of the parameter KL terms that are most influential to model output variability? In our previous work [22], a derivative-based global sensitivity approach has been established to identify unimportant input parameters, for function-valued quantities of interest such as the pressure field. The approach in [22] guides an efficient input dimension reduction strategy, by identifying the KL terms of the input that contribute most to variability of the output field.

5 CONCLUSIONS

We have studied the input and output dimension reduction of elliptic PDEs, with random field input parameters, via the truncated KL expansion technique. In this study, the covariance function of the stochastic process defining the input parameter field is given, and that of the random output is constructed via Monte Carlo sampling. From numerical experiments with both 1D and 2D elliptic PDEs, we observe that when the correlation length is small, very high-dimensional representation is needed to fully resolve the variations in the input field. However, the elliptic operator is not sensitive to high-order KL terms. As a result, the solution of the elliptic PDE only shows strong dependence to the low-order KL terms of the random input field; moreover, the eigenvalues of the solution covariance operator decay very fast. This enables a low-rank representation of the PDE solution in a low-dimensional input parameter space. We then apply these dimension reduction methods in modeling the biotransport process in tumors with uncertain material properties, and demonstrate that the pressure field can be approximated with a low-dimensional representation even for random permeability fields with small correlation lengths. The efficacy of the low-rank ROMs is verified by their capability to recover the PDFs of the pressures at different locations in the flow field.

We demonstrate in this study that the truncated KL expansion can be an effective approach to reduce the output dimensions of an elliptic PDE. This is important for uncertainty quantification of large flow problems with a huge number of spatial dimensions. Although the truncated KL expansion can also reduce the input dimensions, its effect is not apparent when the correlation length of the covariance function is small. Advanced dimension reduction methods, such as global sensitivity analysis and active subspace, need to be developed to tackle input dimension reduction. One example is our recent work on functional derivative-base global sensitivity analysis [22]. More progress will be reported in our future work.

Acknowledgments

M.L. Yu gratefully acknowledge the faculty startup support from the department of mechanical engineering at the University of Maryland, Baltimore County (UMBC).