This paper presents a meta-modeling-based technique to reduce the computational cost of Monte Carlo methods. A stochastic simulation of groundwater flow is substituted with a graph neural network meta-model. This type of neural network can deal with a non-euclidean structure of the input data, which in our case is a graph representing a random field on an unstructured mesh. In order to find the most suitable meta-model, a comparison of the standard support vector regression with spectral and spatial graph convolutional neural networks is provided. Both the Monte Carlo method and the multilevel Monte Carlo method are extended by the meta-model level. In both cases, up to 50% savings in computational cost are achieved while maintaining the accuracy of Monte Carlo estimates.
Groundwater flow in the vicinity of a future nuclear waste repository is the major motivation for our research. In the first stage, our quantity of interest (QoI) is the total flow through the observed area. However, since not all indispensable properties of the rock environment are known, it is not feasible to determine the desired total flow. For that reason, and given the nature of the rock, the missing properties are modeled as random fields (RFs). The probability density function (PDF) of our QoI is under scrutiny. To approximate the PDF, the maximum entropy method is adopted. The method utilizes so-called generalized statistical moments. In order to estimate these moments, the Monte Carlo methods are employed. These methods consist in repeating a random experiment. Depending on the required accuracy, it might result in thousands of groundwater flow simulations, which can significantly affect the total computational cost.
This fact was behind the idea to substitute simulations
with an approximation model, often called a meta-model.
Furthermore, we use the multilevel Monte Carlo method
(MLMC) instead of the standard Monte Carlo method to
reduce the variance of estimates more effectively. Since
the MLMC uses simulations of different accuracy, we aim
to substitute those with the lowest accuracy, performed
in the largest number. Their approximation shouldn’t be
excessively challenging for meta-modeling. We presume
that this approach could reduce the computational cost of
MLMC estimates. It is also necessary to maintain the
accuracy of the estimates to achieve a good PDF
approximation. Graph convolutional neural networks and the support
vector regression [
First, the groundwater flow problem definition and the maximum entropy method brief introduction are provided in section 2. Then, Monte Carlo methods are presented in section 3. Spektral and spatial graph convolutional neural networks used as meta-models are delineated in section 5. Section 6 provides the link between the MLMC and a meta-model, including the related work. Finally, section 7 is devoted to the results and discussion. 2
A 2D benchmark problem of groundwater flow through a
porous medium is used to test our proposed methods. The
same problem can be found in the paper by Blaheta et al.
[
−K∇ p ⋅ n⃗ = 0, K is the hydraulic conductivity, p is the pore water pressure, −K∇ p is Darcy’s velocity, n⃗ is the unit normal vector. We are interested in the total flow Y through the specified area
1
Y = S0 [−K∇ p ⋅ n⃗](1; y)dy.
The problem is prescribed on unstructured meshes and is
solved by the finite element method using the Flow123d
software [
When considering the hydraulic conductivity as a random
field, the total flow Y is a real random variable. The aim is
(1)
(2)
to determine its probability density function r(Y ), which
is nonnegative on a bounded domain W. The maximum
entropy method [
fr(y)r(y)dy = mr; r = 1; :::; R ( ) is an unknown, {mr}rR=1 are values of genwhere r Y
R teiroanlisze(dsemeo[3m, epn.ts1, 3{5f0r]})r=u1seadrefloirnecaarllcyuilnatdienpgengdenenertafluiznecdmoments and satisfying f1 = 1, fr ∈ CR(W), r = 2; : : : ; R. In
R our case, {fr}r=1 are the Legendre polynomials, R ∈ N is the number of moments.
Provided that the system of equations 3 has a solution,
then it has an infinite number of them. Therefore, the
most apposite solution is the one with the maximum
Shannon entropy (defined in [
( ) with the maximum entropy takes the form of finding the global maximum of the functional
H(r) = − SW r(y) ln(r(y))dy (4) under the constraints prescribed in equations 3.
A numerical solution and general limitations of the
MEM are beyond the scope of this paper. A more detailed
insight is provided in [
Since the MEM utilizes moment values {mr}rR=1 that are also random variables derived from some random input X (which is the hydraulic conductivity K in this study). It is necessary to estimate their expected values. Monte Carlo methods are used for this.
The standard Monte Carlo method (MC) is an approach for estimating expected values of some stochastic simulation variable. The basic idea comes from the law of strong (5) (6) (7) (9) where P is a random variable depending on X . According to the central limit theorem: Consequently, the computational cost of reducing the variance (increasing N) of Yˆ can be very high. To overcome this drawback (for details, see [11, p. 4]), the multilevel Monte Carlo method (MLMC) was formulated. 3.1
In the case of the MLMC [
L E[P] = E[P1] + Q E[Pl − Pl−1], l=2 where Pl−1 ≈ Pl . The estimate of the expected value of P improves from P1 to PL. The unbiased estimate of E[P]: where L is the total number of levels, Nl is a number of simulation samples at level l. Input random data xnl are dependent for each simulation pair n at level l. But they are independent across levels. Since Pl − Pl−1 are also independent across levels, the total estimated variance Vˆ of Pˆ has the following form: numbers. Therefore, in order to obtain an unbiased estimate of the expected value Y = E[P(x)], the MC consists in the arithmetic mean of N independent samples where Vˆ1 is an estimate of P1 variance and Vˆl for l > 1 is an estimate of Pl − Pl−1 variance. The MLMC computational cost:
L C = Q NlCl , (10)
l=1 where Cl is a cost of a single simulation sample at level l measured in, for instance, the number of computational operations, execution time, CPU time, etc. Let C1 denote a cost of P1, Cl is a cost of Pl − Pl−1 for l > 1.
Given the MLMC theoretical properties stated in [11, theorem 1], variance Vl should decrease from l = 1 to l = L, while computational cost Cl should increase from l = 1 to l = L.
Our attention is paid to the latter approach, which has the form of finding the minimum of equation 10 under the constraint
Vt = QL1 l NVˆll . (11) = After some calculus and concerning the use of the MEM, Nl are determined with respect to R moment values Nlr = Á¿ÀÁVˆlr ∑iL=1 ¼
Cl
Vt
VˆirCi , r = 1; :::; R, (12) where Vˆlr is an estimated variance of flr(x) − fl−1(x) for r r-th moment at level l. Finally, Nl = r=m1;a::x:;R Nlr.
The mlmc software library [
A meta-model is a simpler and explicit mathematical
function that approximates complicated functions that can be
both implicit and evaluated by a simulation model,
measurements data, or experiments. Alternative names can
be found in the literature, such as surrogate model,
surface response model, emulator, etc. Meta-modeling is the
process of designing meta-models. The common
metamodels used in the finite element analysis include support
vector regression (SVR), artificial neural network (ANN),
Kriging, also known as Gaussian process [
Artificial neural networks have become very popular for
meta-model design. A regression problem is solved by an
ANN trained by supervised learning. The aim is to
obtain an unknown mapping of input neuron activities to an
output neuron activity. With regard to the nature of our
input data, which is a 2D correlated random field, it would
be suitable to use a convolutional neural network (CNN)
meta-model. Nevertheless, CNNs cannot be applied
directly to data on unstructured meshes (see [
There are few ways how to overcome this difficulty. The
basic solution is to directly apply a 1D CNN to nodal
values of an unstructured mesh considered as a vector.
However, it has limited success [
A graph G = (V; E) is an unweighted undirected graph describing the mesh structure on which a random field is generated. V is a set of vertices representing mesh elements, V = {v1; v2; :::; vS}, S is the number of vertices (=mesh elements). The neighborhood of a vertex v is defined as N(v) = {u ∈ V S(v; u) ∈ E}. Each vertex has a feature representing a random field value at the corresponding mesh element. E is a set of edges. An edge connects adjacent mesh elements. An adjacency matrix A is used to represent G on a computer. A variety of real-world problems can be represented as graphs. Imagine social networks, molecules, or in our case, random fields on unstructured meshes.
Graph neural networks (GNNs) are deep learning-based
models that operate on the graph domain [
A graph signal x ∈ RS (vector of all G vertex values) is transformed by a graph Fourier transform from the spatial domain to the spectral domain. Loosely speaking, the standard Fourier transform is connected to the eigendecomposition of the Laplace operator. In the case of graphs, the Laplace operator is represented by the Laplacian matrix [33, p. 4]. The normalized Laplacian matrix L can also be used as a Laplace operator:
L = I − D−1~2AD−1~2; where A is a graph adjacency matrix, and D is a diagonal matrix of vertex degrees. L is a symmetric realvalued positive semi-definite matrix that can be factorized L = U LU T , where U is a matrix of eigenvectors and L is a diagonal matrix of eigenvalues. Then the graph Fourier transform of x is [34, p. 60]: (14) (15) (16) (17) (18) and its inverse:
F(x) = U T x
F−1(xˆ) = U xˆ: Finally, the graph convolution of x and a filter g ∈ Rk: x ∗G g = F−1(F(x) ⊙ F(g)); where ∗G is the convolutional operator on a graph, and ⊙ represents the element-wise Hadamard product. The signal x can be filtered in the spectral domain by a filter gq in this way (for details, see [13, p. 83]):
x ∗G gq = U gq (L)U T x;
where gq (L) is a polynomial of the eigenvalues of the
Laplacian. The filter represents learnable weights.
(19)
(20)
(21)
ChebNet GCN with a Chebyshev convolutional layer was
first introduced by Defferrard et al. [
K−1 gq = Q qkTk( L˜);
k=0 where qk ∈ Rk is a vector of Chebyshev coefficients, ˜ L = 2L~lmax − I, lmax is the maximum eigenvalue from L. The Chebyshev polynomials are defined recursively by Tk≥2(x) = 2xTk−1(x) − Tk−2(x) with T0(x) = 1, T1(x) = x. Then x filtered by g:
K−1 x ∗G g = U ( Q qkTk(L˜ ))U T x: k=0
K−1 x ∗G g = ( Q qkTk(L˜ ))x;
k=0 The eigendecomposition can be avoided [30, p. 10]: where L˜ = 2L~lmax −I. L˜ is called the rescaled graph Laplacian, the eigenvalues are mapped from [0; lmax] to [−1; 1], the Chebyshev polynomials form an orthogonal basis.
Several K settings were tested. Since K > 1 did not bring
improvement for our problem, K = 1 is used in our study.
For this setting, the ChebNet is very similar to the GCN
proposed by Kipf and Welling [
As for spatial GCNs [
Mk(hv(k−1); h(uk−1); xveu)); (22) where hv(k) represents features of vertex v in a hidden layer k, xveu is an optional feature vector of an edge (v; u). Uk and Mk are the update and message functions with learnable parameters.
GraphSage is an aggregation-based learning model
proposed by Hamilton et al. [
A comparison between spectral GCNs and spatial
GCNs is summarized by Wu et al. [30, p. 13]. They
state that spatial GCNs are usually preferred over spectral
GCNs. However, in our problem 2, we aim to approximate
the solution of the elliptic partial differential equation. To
do that numerically, it is possible to use spectral methods
[
There are several approaches how to use meta-models to
reduce the computational cost of Monte Carlo estimates.
The relatively frequent approach (e.g., [
We propose a different approach based on adding a new coarse meta-model-based level to the original Monte Carlo method (MC or MLMC). The MLMC estimate from equation 8 has now the following form:
1 N1 1 1 N2 2 2 Pˆmeta = N1 nQ=1 P˜1(xn) + N2 nQ=1 P1(xn) − P˜1(xn)+ (24) L where P˜1 denotes a meta-model approximation of P . 1 Since the MLMC assumes Pl−1 ≈ Pl , it is completely valid to employ a meta-model as the coarsest level (which we denote as the meta-level) and the difference between P1 and P˜1 as the subsequent first level. Let MLMCmeta denote the multilevel Monte Carlo method with the meta-level.
In order to meet basic theoretical properties of the MLMC (mentioned in section 3.1) the meta-level computational cost C1 should be lower than C2, and meta-level samples variance V1 should be greater than V2. 6.1
The total MLMC cost (see equation 10) depends on Cl , in our case measured in execution time. As a meta-model brings additional computational costs, parts of the learning process, including their costs, are introduced.
First, the number of training samples Ntr is determined, and numerical simulations are executed. Then, for each training sample, a stored random field is pre-processed to become a meta-model input, let Cpr denote the cost of this operation. The pre-processing can slightly differ for different meta-models. The SVR input is a vector of random field elements, while the GCN input is their graph. Afterward, a meta-model is trained, and the cost of this operation is Cml . Once the meta-model is prepared, it is possible to make predictions at the cost of Cpred . In the case of the meta-level, there is no stored random field because there is no simulation executed. Thus a random field is generated and pre-processed at the cost of Cr f .
Let C2 denote the computational cost of a sample at the first MLMCmeta level:
C2 = C1∗ + (C2tr +C2pr(N2 − Ntr))~N2, (25) where C1∗ is the cost of a simulation sample. C2tr is the cost of a meta-model training procedure:
Ctr
2 = CprNtr +Cml +Cpred Ntr; C2pr represents the cost of a first level meta-model prediction sample:
Cpr
2 = Cpr +Cpred .
Let C1 denote the cost of a meta-level sample that utilizes meta-model trained at the first level:
C1 = Cr f +Cpred .
Sample costs Cl for l > 2 are not affected by the metamodeling. Meta-models training is performed on the cluster; 16 CPUs (Intel Xeon Silver 4114 CPU (2.2GHz)) and 16 GB RAM (DDR4 2400 ECC Reg dual rank) are assigned for that task. Also, MLMC simulations are executed in parallel. Therefore, C measured in execution time is not equal to the real elapsed time. However, the savings in C can significantly affect the total elapsed time, especially for a small Vt , which is necessary to obtain a neat PDF by the MEM. (26) (27) (28) 7
The proposed meta-modeling techniques are investigated in this section. The most suitable one is used by Monte Carlo methods. The models are compared for different mesh sizes and random field parameters.
The cross-validation-like procedure is implemented to compare models. Initially, N = 50000 simulation samples are generated. The number of training samples Ntr is empirically determined based on the properties of the MLMC. Since we are interested in Vt ≤ 1e−5 (for the sake of the MEM capability to approximate a PDF neatly), the number of samples at the coarsest level is at least 2000 for our problem. Thus, the procedure is as follows: 2000 training samples (Ntr = 2000) out of N are chosen for metamodel training. The rest is considered to be test data. Validation data accounts for 20% of Ntr. This procedure is repeated 25 times with independent training sets. Given that we use a random field with the exponential correlation function, the logarithm of the RF values is used as meta-models input to facilitate their training. In all of the following cases, the final Nl is determined based on the geometric sequence of the initial number of samples decreasing across levels from N0 L = 100. The 1 = 2000 to N0 number of moments R = 25 is utilized. 7.1
We have tried dozens of GCN topologies. The most
promising ones consist of one ChebNet/GraphSage layer
with 8/60 output channels and the ReLU activation
function, following by a global summation pool and the output
layer with one neuron with the identity activation function.
Table 1 contains main hyperparameters that enable to
obtain an accurate and computationally efficient meta-model.
Another important hyperparameter is the number of output
channels. The optimal number differs across types of GCN
layers. It forms a vector of hidden representations
corresponding to each vertex. It is possible to think of each
channel as responding to some different set of features, so
different channels could become specialized to recognize
different objects as described by Zhang et al. [
Regarding the SVR, the Gaussian radial basis function
(RBF) kernel and the regularization parameter W = 0:06
are used. An explanation of the role of SVR parameters is
provided in [
Three meta-model techniques are compared: the SVR, the ChebNet GCN, and the GraphSage GCN. Recall that low accurate simulations are supposed to be substituted with meta-models. Hence the emphasis is placed on metamodels trained on small meshes. In particular, meshes of 6, 48, and 546 elements are compared. Random field parameters are by default l = 0:1 and s = 1. For a given mesh size, the models are trained and tested on the same data.
Table 2 and Table 3 show the models final accuracy on training data and test data, respectively. The arithmetic mean and the standard error of 25 calculations of the mean squared error (MSE): D1 ∑iD=1(yi −yimeta)2 are provided, where yi is a correct value and yimeta is a predicted value, D is a number of samples. The relative squared error (RSE) is used to compare the models across cases (see Table 4).
Data in the tables show that all models provide similar results in terms of the train MSE and the test MSE, and also, the RSE values are of the same order of magnitude. Importantly, presented data show no significant outliers that would be highly undesirable for our purposes. Although the results are very similar, it can be seen that the optimal meta-model varies depending on the mesh size. While the GraphSage GCN or the SVR is more advantageous for very small meshes, the use of the ChebNet GCN is most suitable in the case of larger meshes.
Nevertheless, additional experiments show that the meta-model approximation is not sufficiently accurate for meshes of thousands of elements. We face the so-called curse of dimensionality due to the limited number of training samples. An increasing number of samples can overcome this difficulty, but many simulations need to be performed in such a case, and learning cost increases. Thus, to keep Ntr = 2000, we limited ourselves to metamodels based on simulations on meshes with a maximum of ca. 1000 elements.
In practice, we see that the ChebNet GCN is not only better for larger meshes (from ca. 500 up to ca. 1000 elements) but also provides a more stable learning process in terms of a smooth decrease of a validation loss than the GraphSage GCN. Therefore the ChebNet GCN is preferred over the GraphSage GCN for our task. 7.3
Random field parameters affect the meta-model learning ability. Changing the standard deviation s only scales features of graph vertices, which our learning procedure can handle. The correlation length l plays a more significant role. As l decreases, the correlation between the features of vertices decreases as well. For small correlation lengths, e.g., l = 0:001, all features are almost uncorrelated, which has a similar effect as increasing the number of vertices with the original correlation length l = 0:1. Thus, increasing the number of training samples is necessary to obtain the same results for l = 0:001.
This reveals one of the SVR drawbacks, which is the requirement of a lot of training samples [2, p. 76]. While the number of training samples 5000 is enough for the ChebNet GCN on meshes of 546 elements and l = 0:001. In the case of the SVR, even 15000 training samples is not enough to get similar results for l = 0:001 as for l = 0:1. For this reason, the ChebNet GCN is preferred for further analysis. It is also worth noting that naturally, fewer training samples are enough for l > 0:1. 7.4
Let use the trained meta-models with the Monte Carlo method described in section 6. To illustrate our approach, we use three different MC with simulations on meshes of 6, 48, and 546 elements, each extended to two-level MLMCmeta. Table 5 shows the ratio between the total cost Cmeta of our MLMCmeta and the total cost C of the original MC, Vt = 1e−5. It is important to note that the ratio Cmeta~C is a bit smaller for Vt << 1e−5, where the metamodel learning cost is negligible due to the larger N2. It emerges that the computational cost savings of at least 50% are achieved for R = 2.
In the rest of this section, MLMCmeta employs the ChebNet GCN on 546 mesh elements. Since the MEM utilizes moment values m , it is important to verify if their estimates by the MLMCmeta are the same as estimates by the original MC. We construct a reference MC (MCref) of 50000 samples on meshes of 217208 elements, moments estimates are denoted as mˆ ire f . Figure 3 shows the MSE: 215 ∑i2=51(mˆ ire f − mˆ i)2 for estimated moment values by the MC, and the MSE: 215 ∑i2=51(mˆ ire f − mˆ imeta)2 for moments estimated by the MLMCmeta.
MC MLMCmeta
Since the error of moments estimates is similar for both the MC and the MLMCmeta, it is a good prerequisite for a fine PDF approximation. Figure 4 shows an example of a PDF approximated by the MEM based on moments from the MC and the MLMCmeta. Given the 25 repetitions, the average KL divergence KL(rre f SSrMC) = 0:034 and KL(rre f SSrMLMCmeta ) = 0:031, for R = 25, Vt = 1e−5. Thus, it is possible to get comparable results by both approaches. Since we cannot effectively train the meta-models on meshes with more than ca. 1000 elements, it is advisable to use initially the MLMC that generally reduces the computational cost compared to the MC. Applying formula 24, a meta-model is trained on the coarsest level, where a mesh should have a small number of elements. Thus, meta-model learning is feasible. To illustrate the MLMC extended by the meta-level, let assume 3 level MLMC with simulations on meshes of 546, 6772 and 87794 elements, and the ChebNet GCN meta-model. RF parameters: l = 0:1, s = 1.
Figure 5 shows the variance decrease from the coarsest level to the finest level. Levels are here defined by a simulation step h. The smaller the h, the finer the simulation. It can be observed that in our current MLMC implementation, the decrease is steeper for lower moments. In general, the more the variance across levels decreases, the more effective the use of the MLMC. Thus, in our case, a higher R results in a less effective MLMC compared to a lower R. The number of moments affects the final Nl and consequently the total computational cost C.
Table 6 provides the ratio between the MLMCmeta total cost Cmeta and the MLMC total cost C. To measure a quality of a PDF approximation, KL divergence KL(rr2e5f SSrMR LMCmeta ) is added to the table, where rr2e5f denotes a PDF based on 25 moments estimated by the MCref, rMR LMCmeta is a PDF approximated from moments estimated by the MLMCmeta. It can be seen that the computational cost savings are greatest for the smallest number of moments. However, R = 2 is insufficient to obtain a decent PDF by the MEM. In our case, we need at least R = 5. On the other hand, with Vt >= 1e−5 and R > 30, the moments estimation error causes the PDF to contain obvious ripples. Using 5 ≤ R ≤ 30, we still achieve savings in the computational cost of at least 10%. For you to get an idea, for R = 25, Vt = 1e−6, the absolute computational costs in seconds are as follows: Cmeta ≈ 73266 and C ≈ 84603.
The effect of R is also manifested in the case of the MC extended by the meta-level. In both cases, it can be noted that cost savings are almost constant for R ≥ 15. 20 15 10 5 r R In this paper, we dealt with the meta-model design for the groundwater flow problem. The motivation was to incorporate a meta-model into the multilevel Monte Carlo method and achieve additional savings in the MLMC computational cost. After comparing the meta-models based on the ChebNet GCN, the GraphSage GCN, and the support vector regression, the ChebNet GCN was selected as the most suitable for our task. The proposed metamodeling technique is effective for simulations on unstructured meshes with a maximum of ca. 1000 elements. Computational cost savings of up to 50% were achieved for both the MC and the MLMC extended by the metalevel. Due to the observed uneven decrease in variances of MLMC estimates across levels, the amount of computational cost savings depends on the number R of generalized statistical moments. In order to obtain a good PDF approximation by the MEM, we required R ≥ 5. In these cases, we were still able to achieve at least 10% savings in computational costs.
Although the obtained results are auspicious, to provide more general conclusions, it is necessary to try our approach with a more complex simulation, which will be more challenging for the meta-model design.
This work was partly supported by the Student Grant Scheme at the Technical University of Liberec through project nr. SGS-2020-3068.
The project leading to this work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 847593.