Learning adaptive coarse spaces of BDDC algorithms for stochastic elliptic problems with oscillatory and high contrast coefficients

In this paper, we consider the balancing domain decomposition by constraints (BDDC) algorithm with adaptive coarse spaces for a class of stochastic elliptic problems. The key ingredient in the construction of the coarse space is the solutions of local spectral problems, which depend on the coefficient of the PDE. This poses a significant challenge for stochastic coefficients as it is computationally expensive to solve the local spectral problems for every realisation of the coefficient. To tackle this computational burden, we propose a machine learning approach. Our method is based on the use of a deep neural network (DNN) to approximate the relation between the stochastic coefficients and the coarse spaces. For the input of the DNN, we apply the Karhunen-Lo\`eve expansion and use the first few dominant terms in the expansion. The output of the DNN is the resulting coarse space, which is then applied with the standard adaptive BDDC algorithm. We will present some numerical results with oscillatory and high contrast coefficients to show the efficiency and robustness of the proposed scheme.


Introduction
Coefficient functions are often one of the main difficulties in modeling a real life problem. Due to the nature of the phenomena, these coefficients can be oscillatory and with high contrast features, all such conditions hinder the development of a fast and accurate numerical solution. Under some conditions, numerical methods like discontinuous Galerkin methods could give robust scheme as in [6,7]. However, more importantly, coefficients are not constant and they usually involve uncertainties and randomness. Therefore, a stochastic partial differential equation (SPDE) is considered for the modeling of these stochastic coefficient functions. There are several popular approaches for solving SPDE numerically, for example, Monte Carlo simulation, stochastic Galerkin method, and stochastic collocation method in [2,3,11]. All these methods require a high computational power in realistic simulations and cannot be generalized to similar situations. In the past decades, machine learning technique has been applied on many aspects such as image classification, fluid dynamics and solving differential equations in [4,17,19,31,12,25,27,5,29]. In this paper, instead of directly solving the SPDE by machine learning, we choose a moderate step that involves the advantage of machine learning and also the accuracy from traditional numerical solver.
We first apply the Karhunen-Loève expansion on the stochastic coefficient as in [10,28,30] to capture the characteristic with finite terms. After this expansion, the SPDE is reduced to a deterministic differential equation and can be solved with a balancing domain decomposition by constraints (BDDC) algorithm, which can prevent the computational problem caused by the refinement of spatial mesh resolution and give an accurate solution rather than just using machine learning technique without suitably designed conditions. After obtaining sufficient number of random solutions from the BDDC method, we can derive the statistics of the original SPDE solution. However, the cost for forming and solving generalized eigenvalue problems in the BDDC algorithm with good accuracy is considerable especially for high dimensional problems. We note that such eigenvalue problems are considered to enrich the coarse space adaptivity in the BDDC algorithm. Therefore, we introduce neural network with a user-defined number of layers and neurons to compensate the increase in computational cost with unexpected times of BDDC solutions.
The adaptive BDDC algorithm with enriched primal unknowns is considered in this paper as its ability on oscillatory and high contrast coefficients has been shown in [13,14]. Among different domain decomposition methods, the considered adaptive BDDC does not require a strong assumption on coefficient functions and the subdomain partitions to achieve a good performance like the standard BDDC algorithm [8,21]. It is because the additional coarse space basis functions computed by the dominant eigenfunctions are related to the ill-conditioning in the standard BDDC algorithm with highly varying coefficients. The introduction of these dominant eigenfunctions thus could greatly improve the robustness of numerical scheme for problems with rough and high contrast coefficients. Moreover, an estimate of condition number of this adaptive BDDC algorithm can be shown to be only controlled by a given tolerance without any extra assumptions on the coefficients or subdomain partitions. For simplicity, in the remaining of the paper, we will call the resulting new numerical scheme combining the machine learning technique and adaptive BDDC algorithm as learning adaptive BDDC algorithm.
Finally, to end this section, we state the main idea of our proposed scheme and the model problem. In this paper, we integrate an artificial neural network with learning abilities into a BDDC algorithm with adaptively enriched coarse spaces for efficient solutions of the neural network approximation of stochastic elliptic problems. Let D be a bounded domain in R d , d = 2, and Ω be a stochastic probability space: where x ∈ D, ω ∈ Ω and ρ(x, ω) is uniformly positive and is highly heterogeneous with very high contrast. Here, we only assume the two-dimensional spatial domain, however, the three-dimensional case is also applicable and is shown to be robust for deterministic elliptic problems in [13]. Moreover, as the model problem (1) can be extended to fluid flow problems, we seldom call the coefficient function ρ(x, ω) as permeability function. The rest of the paper is organized as follows. In Section 2, a brief formulation of the BDDC algorithm with adaptively enriched coarse problems is presented for two-dimensional elliptic problems. Then, the Karhunen-Loève expansion is introduced and the details of the artificial neural network used are clarified in Section 3. Two concrete examples with an explicit analytical expression of the Karhunen-Loève expansion, some network training and testing parameters are presented in the first half of Section 4 and results of various numerical experiments are presented subsequently. Finally, a concluding remark is given in Section 5.

Adaptive BDDC algorithm
The main feature of the adaptive BDDC is the local generalized eigenvalue problem defined on every subdomain interface, which introduces the adaptively enriched primal unknowns. In this section, we repeat the formulation presented by the first and second authors in [13] and give a brief overview of the adaptive BDDC algorithm. We refer to [13,14] for further details, where the analysis of the condition number estimation and the robustness of numerical scheme with oscillatory and high contrast coefficients can be found.

Local linear system
We first introduce a discrete form of the model problem (1) in a deterministic fashion. Let V h be the space of conforming linear finite element functions with respect to a given mesh on D with the mesh size h > 0 and with the zero value on ∂D. We will then find the approximate solution where We assume that the spatial domain D is partitioned into a set of N non-overlapping subdomains {D i }, i = 1, 2, · · · , N , so that D = ∪ N i=1 D i . We note that the subdomain boundaries do not cut triangles equipped for V h . We allow the coefficient ρ(x) to have high contrast jumps and oscillations across subdomains and on subdomain interfaces. Let a i (u, v) be the bilinear form of the model elliptic problem (2) restricted to each subdomain D i defined as In the BDDC algorithm, the original problem (2) is reduced to a subdomain interface problem and solved by an iterative method combined with a preconditioner. The interface problem can be obtained by solving a Dirichlet problem in each subdomain. After choosing dual and primal unknowns on the subdomain interface unknowns, the interface problem is then solved by utilizing local problems and one global coarse problem corresponding to the chosen sets of dual and primal unknowns, respectively. At each iteration, the residuals are multiplied by certain scaling factors to balance the errors across the subdomain interface regarding to the energy of each subdomain problem. The coarse problem aims to correct the global part of the error in each iteration and thus the choice of primal unknowns is important in obtaining a good performance as the number of subdomains increases. The basis for primal unknowns is obtained by the minimum energy extension for a given constraint at the location of primal unknowns and such a basis provides a robust coarse problem with a good energy estimate. We refer to [8,18,21,24] for general introductions to BDDC algorithm.

Notation and preliminary results
To facilitate our discussion, we first introduce some notation. Let S (i) be the Schur complement matrix obtained from the local stiffness matrix A (i) after eliminating unknowns interior to D i , where In the following we will use the same symbol to represent a finite element function and its corresponding coefficient vector in order to simplify the notations.
Recall that X i is the restriction of the finite element space V h to each subdomain D i . Let W i be the restriction of X i to ∂D i . We then introduce the product spaces where we remark that the functions in X and W are totally decoupled across the subdomain interfaces. In addition, we introduce partially coupled subspaces X, W , and fully coupled subspaces X, W , where some primal unknowns are strongly coupled for functions in X or W , while the functions in X, W are fully coupled across the subdomain interfaces.
Next, we present basic description of the BDDC algorithm; see [8,18,21,24]. For simplicity, the two-dimensional case will be considered. After eliminating unknowns interior to each subdomain, the Schur complement matrices S (i) are obtained from A (i) and they form the algebraic problem considered in the BDDC algorithm, which is to find w ∈ W such that where R i : W → W i is the restriction operator into ∂D i , and g i ∈ W i depends on the source term f . The BDDC preconditioner is built based on the partially coupled space W . Let R i : W → W i be the restriction into ∂D i and let S be the partially coupled matrix defined by For the space W , we can express it as the product of the two spaces where W Π consists of vectors of the primal unknowns and W ∆ consists of vectors of dual unknowns, which are strongly coupled at the primal unknowns and decoupled at the remaining interface unknowns, respectively. We define R : W → W such that where R ∆ is the mapping from W to W ∆ and R Π is the restriction from W to W Π . We note that R ∆ is obtained as ∆ is the space of dual unknowns of D i . The BDDC preconditioner is then given by where D is a scaling matrix of the form Here the matrices D i are defined for unknowns in W i and they are introduced to make the preconditioner robust to the heterogeneity in ρ(x) across the subdomain interface. In more detail, D i consists of blocks D where F denotes an equivalence class shared by two subdomains, i.e., D i and its neighboring subdomain D j , and V denotes the end points of F , respectively. We call such equivalence classes F and V as edge and vertex in two dimensions, respectively. In our BDDC algorithm, unknowns at subdomain vertices are included to the set of primal unknowns and adaptively selected primal constraints are later included to the set after a change of basis formulation. For a given edge F in two dimensions, the matrices D F are chosen simply as multiplicity scalings, i.e., 1/2, but for a general case when ρ(x) has discontinuities, different choice of scalings, such as ρ-scalings or deluxe scalings, can be more effective. The scaling matrices D (l) V can be chosen using similar ideas. We refer to [16] and references therein for scaling matrices.
We note that by using the definitions of R and S, the matrix in the left hand side of (4) can be written as In the BDDC algorithm, the system in (4) is solved by an iterative method with the preconditioner (5). Thus its performance is analyzed by estimating the condition number of

Generalized eigenvalue problems
After the preliminaries, we are ready to state the generalized eigenvalue problems which introduce the adaptive enrichment of coarse components. For an equivalence class F shared by two subdomains D i and D j , the following generalized eigenvalue problem is proposed in [15]: where S (i) F are the block matrices of S (i) and D i corresponding to unknowns interior to F , respectively. The matrix S (i) F are the Schur complement of S (i) obtained after eliminating unknowns except those interior to F . In addition, for symmetric and semi-positive definite matrices A and B, their parallel sum A : B is defined by, see [1], where (A + B) + is a pseudo inverse of A + B. We note that the problem in (7) is identical to that considered in [9] when D (i) F are chosen as the deluxe scalings, i.e., We note that a similar generalized eigenvalue problem is considered and extended to three-dimensional problems in the first and second authors' work [13].

Learning adaptive BDDC algorithm
In addition to the BDDC algorithm, we perform the Karhunen-Loève (KL) expansion to decompose the stochastic permeability ρ(x, ω) for the preparation of the learning adaptive BDDC algorithm. As shown in [30], only a few terms in the KL expansion are already enough to approximate the stochastic permeability with a reasonable accuracy, which can be efficiently reduce computational cost. Once we obtain the KL series expansion, we can easily produce a certain amount of sample data as a training set for the artificial neural network. After the training process, the neural network captures the main characteristics of the dominant eigenvectors from the adaptive BDDC algorithm, and neural network approximation can be obtained. In our proposed algorithm, we still need some resulting samples from the deterministic BDDC algorithm to train the neural network, which is a process of a relatively high computational cost. However, after this setup, we can apply this network to obtain an approximate solution directly without the use of other stochastic sampling methods such as Monte Carlo simulation, which involves even higher computational cost for the generation of large number of random samples or realizations.
In this section, we discuss the proposed algorithm in details. Moreover, our resulting neural network can be also applied on problems with similar stochastic properties. It thus saves the cost of training again and we can readily use the trained neural network for prediction. Some testing samples will be generated to test the network performance and to show its generalization capacity in the later section.

Karhunen-Loève expansion
To start our learning adaptive algorithm, we first apply the KL expansion on the stochastic field ρ(x, ω). To ensure the positivity of stochastic permeability field, a logarithmic transformation is considered as K(x, ω) = log(ρ(x, ω)). Let C K (x,x) be the covariance function of K(x, ω) at two locations x andx. Since the covariance function C K (x,x) is symmetric and positive definite, it can be decomposed into: where λ i and f i are eigenvalues and eigenfunctions computed from the following Fredholm integral equation: where {f i (x)} are orthogonal and deterministic functions. It is noted that for some specially chosen covariance functions, we can find the (λ, f ) eigenpair analytically. We will list some examples in the section of numerical results. After solving the Fredholm integral equation, we are ready to express the KL expansion of the log permeability: Here E[K](x) is the expected value of K(x, ω) and ξ i are identically independent distributed Gaussian random variables with mean 0 and variance 1. As mentioned above, a few terms in KL expansion can give a reasonably accurate approximation, which can be explained now as we can just pick the dominant eigenfunctions in (11) and the other eigenvalues will decay rapidly. Therefore, KL expansion is an efficient method to represent the stochastic coefficient and flavour the subsequent process. Although analytical eigenfunctions and eigenvalues cannot always be found, there are various numerical methods to compute the KL coefficients in [23,26]. After the generation of some sets of {ξ i }, which are the only stochastic variables in the permeability after KL expansion, as the input of neural network, we plug back KL expanded permeability coefficients into the adaptive BDDC algorithm, and the resulting dominant eigenfunctions in the coarse spaces will be the target output of the our designated neural network.

Neural network
There are many types of neural networks for different usages. Since in our learning adaptive BDDC algorithm the desired neural network is to capture numeric feature from the data as a supervised learning, a fully connected feedforward neural network is chosen. An illustration of the neural network structure is shown in Figure 1. In a general feedforward neural network, there are L hidden layers between the input and output layer, and the arrows from layers to layers indicate that all the neurons in the previous layer are used to produce every neuron in the next layer. Computational cost increases with L, but there are no significant improvements in predictions of our tests, L = 1 is thus chosen throughout numerical experiments.
We denote the number of neurons in the input layer to be R, which is the number of terms in the truncated KL expansion. The inputs are actually the gaussian random variables {ξ i } in one realization. The output of the network is a column vector that consists of all dominant eigenvectors in the coarse space. The number of neurons in the output layer is denoted as O, which is the sum of length of all resulting eigenvectors obtained from the deterministic BDDC algorithm. This value O depends on many factors, such as the geometry of the spatial domain, the structure of the grid, and parameters of the BDDC algorithm.
In the training, we use the scaled conjugate gradient algorithm in [22] to do the task of minimization. Therefore, unlike the gradient descent method, we do not need to determine the value of learning rate at each step. However, the initial choice for the scaled conjugate gradient algorithm may still affect the performance of training. Before training, we have to generate training samples and decide some conditions on stopping criteria. When the minimization process is completed for all the training samples, we are said to finish one epoch. There are two stopping criteria that are usually considered. They are the minimum for the cost function gradient and the maximum number of epochs trained. When either of these conditions are satisfied, the training is stopped and the network will be tested for performance on a testing set.
In both training and testing results, in order to obtain an error measure less affected by scales of data set, a normalized root mean squared error (NRMSE) is used as an error of reference to see (2) . . . . . .   how well are the network estimation and prediction: where ξ (i) ∼ N (0, 1) ∈ R R is a column vector of random variables that follow the standard normal distribution in the i-th realization, y(·) is the dominant eigenvectors obtained from the adaptive BDDC algorithm when the KL expanded stochastic permeability function is used, and f N N (·) is the regression function describing the neural network (NN). Here, O is the dimension of network output, and M is the number of samples. For simplicity, the considered NRMSE can be understood as the root mean squared error divided by the maximum element of the absolute difference of all the eigenvectors. Overall, the proposed scheme can be concluded as follows: • Step 1: Perform KL expansion on the logarithmic stochastic permeability function K(x, ω)

Numerical results
In this section, supporting numerical results are presented to show the performance of our proposed learning adaptive BDDC algorithm. We will consider various choices of permeability coefficients ρ(x, ω) with two major stochastic behaviours. As the concern of our numerical tests are on the learning algorithm, we fix all the parameters that are not related. In all the experiments, D ∈ R 2 is chosen to be a unit square spatial domain, and it is partitioned into 16 uniform square subdomains with H = 4 −1 as the coarse grid size. Each subdomain is then further divided into uniform grids with a fine grid size h = 32 −1 . In the following, we will first describe the considered stochastic coefficients and training conditions.

Choices of stochastic coefficients
From the KL expansion expression (11), we know the eigenvalues and eigenfunctions computed from the Fredholm integral equation are related to the covariance function of K(x, ω), which can be treated as the starting point of KL expansion. Throughout the experiments, two specially chosen covariance functions whose eigenvalues and eigenfunctions in the Fredholm integral equation can be computed analytically are considered with different expected values. For x = (x 1 , x 2 ) T and x = (x 1 ,x 2 ) T , these two covariance functions are 1. Brownian sheet covariance function: where the corresponding eigenvalues and eigenfunctions are 2. Exponential covariance function: where σ 2 K and η i are the variance and correlation length in the x i direction of the process, respectively. The corresponding eigenvalues and eigenfunctions are where r 1,i is the i-th positive root of the characteristic equation (r 2 1,i η 2 1 − 1) sin(r 1,i ) = 2η 1 r 1,i cos(r 1,i ) in the x 1 direction. And the same for r 2,j , which is in the j-th positive root of the characteristic equation in the x 2 direction. After arranging the roots r 1,i and r 2,j in ascending order, we can immediately obtain a monotonically decreasing series of λ k , which flavours our selection of dominant eigenfunctions in the truncated KL expansion with R terms.
For clarification, the positive indices i and j are mapped to index k so that λ k is monotonically decreasing. It is noted these eigenvalues and eigenfunctions are just the product of solutions of the one-dimensional process in the Fredholm integral equation, since the covariance functions are separable.

Training conditions
In the training, there are different sets of parameters that could affect the performance of the neural network. We list down as follows: After we obtain an accurate neural network from the proposed algorithm, besides the NRMSE of testing set samples, we also present some characteristics of the preconditioner based on the approximate eigenvectors from the learning adaptive BDDC algorithm, which include the number of iterations required, minimum and maximum eigenvalues of the preconditioner. To show the performance of preconditioning on our predicted dominant eigenvectors, we will use the infinity norm to measure the largest difference and also the following symmetric mean absolute percentage error (sMAPE), proposed in [20], to show a relative error type measure with an intuitive range between 0% and 100%: where q i are some target quantities from preconditioning and q i is the same type of quantity as q i from preconditioning using our predicted eigenvectors. In the following subsections, all the computation and results were obtained using MATLAB R2019a with Intel Xeon Gold 6130 CPU and GeForce GTX 1080 Ti GPU in parallel.

Brownian sheet covariance function
The first set of numerical tests bases on the KL expansion of Brownian sheet covariance function with the following expected functions E[K](x): (A) 10 s (B) 5 sin(2πx 1 ) sin(2πx 2 ) + 5.
Expected function A is a random coefficient as s ∈ [−1, 1] is randomly chosen for each fine grid element, and expected function B is a smooth trigonometric function with values between [0, 10]. The common property between these two choices is that the resulting permeability function ρ(x, ω) is highly oscillatory with magnitude order about 10 −2 to 10 5 , which are good candidates to test our method on oscillatory and high contrast coefficients. Here, we show the appearance of mean functions A and B and the corresponding permeability function K(x, ω) in Figure 2. We can observe that from the first row of Figure 2 to the second row, that is, from the appearances of expected function E[K](x) to the logarithmic permeability coefficient K(x, ω), there are no significant differences except only a minor change on the values on each fine grid element. Nevertheless, these small stochastic changes cause a high contrast ρ(x, ω) after taking exponential function as in the third row of Figure 2.
In Table 1, the training results of these two expected functions A and B under the training conditions mentioned are presented. Although the random initial choice of conjugate gradient algorithm for minimization may have influences, in general, the neural network of expected function A usually needs more epochs until stopping criterion is reached. However, this increase in training epochs does not bring a better training NRMSE when compared with expected function B. We can see the same phenomenon when the testing set is considered. The main reason is due to the smoothness of expected function B, the neural network can capture its characteristics easier than a random coefficient.

Case
Epochs trained Training NRMSE A 2.29e+05 1.97e-02 B 5.39e+04 3.77e-03 In the testing results, we consider another set of data containing M = 500 samples as a testing set. After obtaining the predicted eigenvectors from neural network, we plug it into the BDDC preconditioned solver and obtain an estimated preconditioner. For the estimated preconditioner, we report several properties such as the number of iteration needed for the iterative solver, the maximum and minimum eigenvalues of the preconditioned system. Therefore, to better show the performance of network, besides the testing error NRMSE, we also use the sMAPE and l ∞ norm of difference of the iteration numbers, the maximum and minimum eigenvalues of the estimated preconditioner and the target preconditioner. The comparison results are listed below:

Case
Testing NRMSE sMAPE ( l ∞ error) in Iteration number λ min λ max A 6.58e-02 6.94e-02 ( 2 ) 2.58e-07 ( 4.02e-06 ) 4.45e-03 ( 3.58e-02 ) B 6.73e-03 6.93e-02 ( 1 ) 3.19e-05 ( 1.50e-04 ) 2.01e-03 ( 4.94e-02 ) It is clear that the testing NRMSE of expected function B is also much smaller than expected function A. One possible reason is because the function used is smooth, and the neural network can better capture the property of the resulting stochastic permeability function. We can see that values of the testing NRMSE of both expected function A and expected function B are just a bit larger than the training NRMSE. This suggests that the neural network considered is capable to give a good prediction on dominant eigenvectors in the coarse space even when the magnitude of coefficient function value changes dramatically across each fine grid element. Nevertheless, even though a more accurate set of eigenvectors can be obtained for the expected function B, but its performance after the BDDC preconditioning is not always better than that of expected function A. Obviously, l ∞ errors in the aspects of λ max and λ min are larger for the expected function B.
Detailed comparisons can be found in Figure 3, where the area of overlapping represents how close is the specified quantity of the estimated preconditioner compared to the target preconditioner. In the first row of figures, we can see the number of iteration required for the predicted eigenvectors is usually more than the target iteration number, however, the differences in the iteration number is not highly related to the difference in minimum and maximum eigenvalues. We can see more examples to verify it in this and later results. The bin width of the histogram of minimum eigenvalues for expected function A is about 10 −7 , which is in coherence with the l ∞ error, but we eliminate this small magnitude effect and the sMAPE is considered. The high ratio of overlapping area for expected function A confirms the usage of sMAPE as a good measure, as the resulting sMAPE 2.57e-07 of expected function A is much smaller than 3.19e-05 of expected function B.
For the last column of errors, the l ∞ norm of maximum eigenvalues is larger for expected function B. Nevertheless, it does not mean that the neural network of expected function B performs worse. From its histogram, we can see an extremely concentrated and overlapping area at the bin around 1.04, moreover, there are some outliers for the prediction results, which is one of the source for a larger l ∞ norm. Therefore, the sMAPE can well represent the performance of preconditioning on the estimated results, and we can conclude that both neural networks show a good results on oscillatory and high contrast coefficients, which can be represented by NRMSE and sMAPE. In the next subsection, besides an artificial coefficient, we will consider one that is closer to a real life case.
Although the l ∞ norm seems not describe the characteristics of performance in a full picture, it is still a good and intuitive measure, for example, in the difference of iteration number, we can immediately realize how large will be the worst case.

Exponential covariance function
To test our method on realistic highly varying coefficients, we consider expected functions that come from the second model of the 10th SPE Comparative Solution Project (SPE10), with the KL expanded exponential covariance function as the stochastic source. For clarification, in this set of experiment, we first use the modified permeability of Layer 35 in the x-direction of SPE10 data as the expected function E[K](x) with (σ 2 K , η 1 , η 2 ) = (1, 0.25, 0.25) in the exponential covariance function C K (x,x). We then train the neural network and the resulting neural network is used to test its generalization capacity on the following different situations: 1. Different stochastic behaviour: • All remain unchanged except (η 1 , η 2 ) = (0.2, 0.125) 2. Different mean permeability function: • The expected function is changed to the Layer 34 in the x-direction of SPE10 data, but the stochastic parameters are unchanged, i.e., (σ 2 K , η 1 , η 2 ) = (1, 0.25, 0.25) Note, all the permeability fields used are modified into our computational domain and shown in Figure 4 for a clear comparison. In the Layer 35 permeability field, four sharp properties can be observed. There are one blue strip on the very left side, and a reversed but narrower red strip is next to it. In the top right corner, a low permeability area is shown and just below it, a small high permeability area is located in the bottom right corner. So a similar but with slightly different properties of Layer 34 are chosen to test the generalization capacity of our neural network.
Before showing the training and testing records of neural network, we further present the logarithmic permeability K(x, ω) and permeability ρ(x, ω) when (σ 2 K , η 1 , η 2 ) = (1, 0.25, 0.25) and  Figure 5. The sets of {ξ i } used in these two column of realizations are the same. We can thus see the change in parameters does cause the stochastic behaviour to change, too. We present the training and testing records of neural network as below in Table 3 and Table 4, where we specify the expected functions considered in the first column.

Case
Epochs trained Training NRMSE Layer 35 2.62e+05 1.52e-02  In Table 4, we clarify that the results of Layer 35 * also use Layer 35 permeability as the expected function but with (η 1 , η 2 ) = (0.2, 0.125). Each row corresponds to the results of different testing sets, however, all the results are obtained using the same neural network of Layer 35 in Table 3. Therefore, the testing results of Layer 35 are used as a reference to decide whether the results of Layer 35 * and Layer 34 are good or not.
We first focus on the testing results of Layer 35. From the corresponding histograms in the left column of Figure 6, besides the iteration number required for predicted eigenvectors are still usually larger by 1, we can observe a large area of overlapping, in particular, there is only about 10% of data that are not overlapped in the histogram of maximum eigenvalues. Moreover, the small NRMSE and sMAPE again confirm with the graphical results. For generalization tests, we then consider two other testing sets with different stochastic behaviour and expected function. When the stochastic behaviour of Layer 35 is changed to (σ 2 K , η 1 , η 2 ) = (1, 0.2, 0.125), we list the corresponding sMAPE and l ∞ norms of differences in iteration numbers, maximum eigenvalues and minimum eigenvalues of preconditioners in the second row of Table 4. We can see the values of errors are very similar to the first row and the NRMSE is even smaller. Therefore, we expect that the Layer 35 neural network can generalize the coefficient functions with similar stochastic properties, and give a good approximation on dominant eigenvectors for the BDDC preconditioner. It is then verified by the histograms in the right column of Figure 6. Moreover, a large population of samples are concentrated around 0 in the histograms of difference between the estimated preconditioner and target preconditioner in Figure 7. All these show our method performs well with stochastic oscillatory and high contrast coefficients, and also coefficient functions with similar stochastic properties.
Finally, we compare the performances when the expected function is changed to Layer 34. Note, the NRMSE, sMAPE and l ∞ errors of Layer 34 all are much higher than those of Layer 35, which are crucial clues for a worse performance. Although the permeability fields of Layer 34 and Layer 35 share some common properties and have a near magnitude, unlike previous results for Layer 35 and Layer 35 * , we can see in Figure 8, the prediction output and target output on maximum and minimum eigenvalues of preconditioners are even two populations with different centres and only a little overlapping area. The main reason is a minor difference in K(x, ω) can already cause a huge change after taking exponential function. Therefore, when E[K](x) is changed to another entirely different mean permeability with just some common properties, the feedforward neural network may not be able to capture these new characteristics, which will be part of our future research.

Conclusion
A new learning adaptive BDDC algorithm is introduced and it consists of three main parts, which are the Gaussian random variables in Karhunen-Loève expansion, artificial neural network and the dominant eigenvectors obtained from the adaptive BDDC algorithm in the coarse spaces. The considered neural network acts as a fast computational tool, which turns the input Gaussian random variables into predicted dominant eigenvectors as output. In addition, the neural network in the proposed algorithm is suitable for other permeability coefficients with similar stochastic properties for generalization purpose. Numerical results confirm the efficiency of the proposed algorithm and the generalization abilities. On the other hand, we currently just use the simplest feedforward neural network structure without retraining when new characteristics are entered. In our future plan, we will improve the neural network and include deep learning technique for higher accuracy in prediction and preconditioning, and also better generalization ability on wider aspects of applications.