Generalized Neumann Expansion and Its Application in Stochastic Finite Element Methods

An acceleration technique, termed generalized Neumann expansion (GNE), is presented for evaluating the responses of uncertain systems. The GNE method, which solves stochastic linear algebraic equations arising in stochastic finite element analysis, is easy to implement and is of high efficiency. The convergence condition of the new method is studied, and a rigorous error estimator is proposed to evaluate the upper bound of the relative error of a given GNE solution. It is found that the third-order GNE solution is sufficient to achieve a good accuracy even when the variation of the source stochastic field is relatively high. The relationship between the GNEmethod, the perturbationmethod, and the standard Neumann expansionmethod is also discussed. Based on the links between these three methods, quantitative error estimations for the perturbationmethod and the standard Neumannmethod are obtained for the first time in the probability context.


Introduction
Over the past few decades, the stochastic finite element method (SFEM) has attracted increasing attention from both academia and industries.The objective of SFEM is to analyze uncertainty propagation in practical engineering systems that are affected by various random factors, for example, heterogeneous materials and seismic loading.The current research of SFEM has two main focuses: probabilistic modeling of random media and efficient numerical solution of stochastic equations.Properties of heterogeneous materials (e.g., rocks and composites) vary randomly through the medium and are generally modeled with stochastic fields.Various probabilistic methods [1][2][3][4][5][6][7] have been developed to quantitatively describe the random material properties, which include the spectral representation method [1,2], the Karhunen-Loève expansion method [3,4], and the Fourier-Karhunen-Loève expansion method [5,6], to name a few.In these methods, the random material properties are represented by a linear expansion such that they can be directly inserted into the governing equations of the physical system, which facilitates further numerical analysis.The resulting equations are stochastic equations, whose solution is the second main focus of SFEM research and also the focus of this work.
For static and steady-state problems, the resulting SFEM equations have the following form [8][9][10][11][12][13]: where A 1 , . . ., A  are  ×  symmetric and deterministic matrices, u() the unknown stochastic nodal vector, and f() the nodal force vector.If the material properties are modeled as Gaussian fields, then  1 ,  2 , . . .,   are independent Gaussian random variables.The uncertainties from the right-hand side of (1) are relatively easy to handle because they are not directly coupled with the random matrix from the left-hand side.
A number of methods have been developed for the solution of (1), including Monte Carlo methods, the perturbation method, the Neumann expansion method, and various projection methods.The Monte Carlo methods [14][15][16][17][18][19] first generate a number of samples of the random material properties, then solve each sample with the deterministic FEM solver, and finally evaluate the statistics from these deterministic solutions.The Monte Carlo methods are simple, flexible, friendly to parallel computing, and not restricted by the number of random variables [20].The computational cost of Monte Carlo methods is generally high especially for large-scale problems, although there exist a number of efficient sampling techniques such as important sampling [14,19], subset simulation [15,17], and line sampling [18,19].Owning to its generosity and accuracy when sufficient sampling can be ensured, the Monte Carlo method is often taken as the reference method to verify other SFEM techniques.The perturbation method [11,[21][22][23][24][25] expresses the stiffness matrix, the load vector, and the unknown nodal displacement vector in terms of Taylor expansion with respect to the input random variables, and the Taylor expansion coefficients corresponding to the unknown nodal displacement are then solved by equating terms of the same order from both sides of (1).Its accuracy decreases when the random fluctuation increases, and there has not been a quantitative conclusion for the relationship between the solution error and the level of perturbation, making it difficult to confidently deploy the perturbation method in practical problems.The Neumann expansion method [9,26,27] is an iterative process for the solution of linear algebraic equations.Because the spectral radius of the iterative matrix must be smaller than one, it is often considered not suitable for problems with large random variations.However, the fluctuation range allowed by the Neumann expansion method is generally higher than the perturbation methods, because it is easier to adopt higher order expansions in the former one.It should be noted that the stationary iterative procedure employed in Neumann expansion is not as efficient as the preconditioned conjugate gradient procedure employed in Monte Carlo methods.Therefore, the computational efficiency of the Neumann expansion method is inferior to direct Monte Carlo simulations.In the projection methods [3,4,12,[28][29][30][31][32][33][34][35], the unknown nodal displacements are projected onto a set of orthogonal polynomials or other orthogonal bases, after which the Galerkin approach is employed to solve the projection coefficients.However, the projection methods typically result in an equation system much larger than the deterministic counterpart, and thus these methods are often limited to cases with a small set of base functions and a few independent random variables.In order to handle problems with more random variables, the stochastic collocation methods [36][37][38][39][40][41][42] have been developed, in which a reduced random basis is used in the projection operation and the projection coefficients are determined by approximating the precomputed Monte Carlo solutions.However, the stochastic collocation methods are believed to be less rigorous and their accuracy is poorer than the full projection methods using the Galerkin approach.Other variants of the projection procedures can be found in [13,43,44].The fatal drawback of the projection types of solvers is their low efficiency [45], because the size of the base function set grows exponentially with the number of random variables and the order of the polynomials used.Moreover, the accuracy of the projection methods is criticized by [46,47].For linear stochastic problems, a joint diagonalization solution strategy was proposed in [8][9][10], in which the stochastic solution is obtained in a closed form by jointly diagonalizing the matrix family in (1).This method is not sensitive to the number of random variables in the system, while there has not been a mathematically rigorous proof for its accuracy.
We propose a generalized Neumann expansion method for the solution of (1).The new GNE method is many times more efficient than the standard Neumann method, which is compromised by increased memory demand.In the context of probability, a mathematically rigorous error estimation is established for the new GNE scheme.It is also proved that for (1), the GNE scheme, the standard Neumann expansion method, and the perturbation method are equivalent.Hence, the stochastic error estimation obtained here also shed a light on the standard Neumann method and the perturbation method, which has remained as an outstanding challenge in the literature.The rest of the paper is organized as follows.The GNE scheme and its error estimation are derived in Section 2, after which it is validated and examined in Section 3 via a comprehensive numerical test.Concluding remarks are given in Section 4.

Generalized Neumann Expansion Method
Consider a linear elastic statics system with Young's modulus modeled as a random field.By using the Karhunen-Loève expansion [3,4] or the Fourier-Karhunen-Loève expansion [5,6], random Young's modulus can be expressed as where (x, ) denotes random Young's modulus at point x,  0 (x) is the expectation of (x, ),   () are uncorrelated random variables, and   and   (x) are deterministic eigen values and eigen functions corresponding to the covariance function of the source field (x, ).If (x, ) is modeled as a Gaussian field,   () become independent standard Gaussian random variables.For the sake of simplicity, we drop the random symbol  and for the rest of the paper use   to denote the random variables in (2).Following the standard finite element (FE) discretization procedure, a stochastic system of linear algebraic equations can be built from (2) In the above equation, A 0 is the conventional FE stiffness matrix, which is symmetric, positive definite, and sparse while A  are symmetric and sparse matrices.The deterministic matrices A  share the same entry pattern as the stiffness matrix A 0 , but they are usually not positive definite.

Related Methods.
The Monte Carlo methods first generate a set of samples of   in (3), and each sample of   corresponds to a deterministic equation, which can be solved by the standard FEM method.
In the perturbation procedure, the unknown random vector u is expanded as [21,25,26] Substituting ( 4) into (3) yields from which the components u 0 , u  , u  of u are solved as follows: It should be noted that, in the stochastic context, there has not been a rigorous convergence criterion for the perturbation scheme described above.
In the Neumann expansion method, the solution of ( 3) is written as where ΔA = ∑    A  , B = A −1 0 ΔA and I is the identity matrix.According to Neumann expansion solution ( 7) can be simplified as where u 0 = A −1 0 f.The spectral radius of matrix B must be smaller than 1 to ensure the convergence of (9).8) can be generalized to the inverse of multiple matrices:

The Generalized Neumann Expansion Method. The Neumann expansion (
The generalized Neumann expansion (GNE) in ( 10) can be proved as follows: where  is the total number of B  .
The Neumann expansion ( 8) and the generalized Neumann expansion (10) are equivalent in mathematics and only differ in organizing the matrix operations.This difference can lead to a significant acceleration in computing the solution of (3).Specifically, let B  = A −1 0 A  , and applying the GNE in ( 10) to (3) gives in which Equation ( 12) is the GNE solution.

Convergence Analysis and Error Estimation.
The convergence condition of the GNE ( 10) is essentially the same as the Neumann expansion (8), that is, (∑   B  ) < 1, where (⋅) denotes the spectral radius of a matrix.Hence, for the GNE solution in (12), the convergence condition is It should be noted that where ‖ ⋅ ‖  2 denotes the spectral norm of a matrix, and the symmetric and positive definite matrix C is defined as As a type of matrix norm, the spectral norm has the following properties: where P and Q are general matrices and ℎ is a positive constant.Hence, Based on the above relation, a sufficient but not necessary condition for convergence is The maximum spectral radius  max = max((Α −1 0 A  )) can be obtained from A 0 and A  .Equation ( 19) can be used to verify the convergence of GNE in advance.When   are independent standard Gaussian random variables, their absolute values are smaller than 3 with the probability of 99.7%.In this case, (19) can be further simplified: Considering the th-order GNE (10) or Neumann expansion (8), their error can be estimated as follows: where   is the spectral norm of B, which is assumed to satisfy The exact solution of ( 3) is where The th-order GNE or Neumann expansion solution is given by The norm of the solution error is defined as where   is the spectral radius of B.
The relative error of the solution is defined as Comparing ( 21) and ( 26), it can be seen that the accuracy for the solution vector u is even better than the accuracy for the matrix inverse (I + B) −1 .

Solutions for More General Stochastic Algebraic Equations.
For the solution of (3), the GNE solution in (12) shares the same form as the perturbation solution in (5).Hence, for solving (3), the GNE method, the perturbation method, and the standard Neumann expansion method are identical.The error estimation obtained in (26) also holds for the perturbation approach.However, for the following more general stochastic algebraic equations [26]: the three methods are not equivalent.Specifically, applying the perturbation method to the above equation yields Equation ( 28) can be reorganized as If   =   =   , then the perturbation components u 0 , u  , u  , . . ., u  1  2 ⋅⋅⋅  of u can be solved as Solving (27) with the Neumann expansion gives where in which Because higher order entries are included in B, the term BA −1 0 (f 0 + Δf) (the first-order term of the Neumann solution in (32)) does not correspond to the first-order entries of u.Therefore, the Neumann solution in (32) and the perturbation solution in (30) are different for the general stochastic equation (27).
Solving (27) with the GNE approach gives where In summary, when solving (27), the GNE method and the Neumann method are identical, while they are not equivalent to the perturbation procedure.However, if the higher order terms in the GNE or Neumann expansion solutions are omitted, they degenerate into the perturbation solution.Two simplified examples are given here to demonstrate the difference between the GNE and the perturbation method.

Example I.
Only first-order entries are included in both sides of ( 27) such that the equation becomes The first-order perturbation solution is and the first-order GNE solution is Comparing ( 36) and ( 37), it can be seen that the GNE solution contains an extra second-order term − ∑  ∑      A −1 0 A  f  .

Example II.
Up to second-order terms are included in the left-hand side of ( 27), while only the constant term is included in the right-hand side.In this case, the equation becomes Denote u 0 = A −1 0 f 0 , B  = A −1 0 A  , and B  = A −1 0 A  ; then the first-order perturbation solution can be expressed as while the first-order GNE solution is Comparing ( 39) and ( 40), the second-order terms − ∑  ∑      A −1 0 A  u 0 in the GNE solution do not appear in the perturbation solution (39).

Example Description.
As shown in Figure 1, a plain strain problem with an orthohexagonal cross section is considered here.The length of each edge of the hexagon is 0.4 m.The Poisson ratio of the material is set as 0.2, and Young's modulus is assumed to be a Gaussian field with a mean value of 30 GPa and a covariance function defined as in which  1 =  1 −  2 and  2 =  1 −  2 are the coordinate differences and  1 and  2 are adjustable constants.As shown in Figure 2, the cross section is discretized with triangular linear elements, where each side is divided into 8 equal parts.There are in total 384 elements and 217 nodes.The relative error is defined as in which u MC is the nodal displacement vector solved with the Monte Carlo method and u GNE is the GNE solution.Unless stated explicitly, the numerical testing is performed for 4340 samples generated from the Gaussian field defined in (41).
The convergence condition of the GNE method has been rigorously proved.The main purpose of this example is to examine the error distribution and its relationship with various random parameters, for which a number of numerical tests are performed in the following subsections.

Error Estimation.
The third-order GNE solution is adopted.The constants  1 and  2 in (41) are varied to examine the error estimator in (26).The constant  1 is the variance of the Gaussian field, and the constant  2 controls the shape of the covariance function.The effective correlation length  eff of the covariance function () is defined by For any two points in the medium, if the distance between them is larger than  eff , their material properties are only weakly correlated.The effective correlation length of the Gaussian field defined in (41) is In the first set of tests, the constant  1 is varied to change the variance of the Gaussian field.The parameter settings are In the second set of tests, the constant  2 is varied to change the shape of the covariance function.The parameter settings are 3.2.1.The Influence of  1 .Corresponding to  1 = 0.01 and  1 = 0.04, the standard deviation  of the stochastic field is 0.1 and 0.2, respectively.A total number of 4340 samples are generated for each case.For each realization, both the Monte Carlo solution and the third-order GNE solution are computed.The relationship between the spectral radius and the relative error is drawn in Figures 3 and 4, in which the red crosses show the actual relative errors corresponding to each sample and the blue line shows the theoretical error estimation.It can been seen that for both cases, the theoretical error estimator in (26) holds accurately and its effectiveness is not affected by the change of the variance controlled by  1 .

The
Influence of  2 .The constant  2 is adjusted to take three different values 0.8 × 10 −2 , 2.4 × 10 −2 , and 4.0 × 10 −2 , while the constant  1 remains unchanged.Following (41), the effective correlation length grows with the increase of  2 .Again, 4340 samples are generated for each stochastic field, and for each realization both the Monte Carlo solution and the third-order GNE solution are computed.Figure 5 plots the relationship between spectral radius and the relative error, which verifies that, for all cases, the theoretical error estimation in (26) holds.The changing of  2 , which changes the variation frequency of the stochastic field, does not affect the error estimation.
As indicated in (44), the effective correlation length  eff is directly associated with the parameter  2 .To investigate the influence of  eff on the relative error, the empirical cumulative distribution functions (ECDF) of relative errors are plotted in Figure 6 for different  eff cases.Comparing the three ECDF curves at different effective correlation lengths, it is found that with the decrease of the effective correlation length, the accuracy of the GNE solution improves.

The Relative Error for Other Types of Stochastic Fields.
In the last two subsections, the error estimation analysis has taken the Gaussian assumption; that is, all random variables   in (3) are assumed to be independent Gaussian random variables.In this subsection, we investigate the relative error of the GNE method for other types of stochastic fields.Specifically, the random variables   in (3) are assumed to have the two-point distribution, the binomial distribution, and the  distribution, respectively.For the case of two-point distribution, values {−0.3 0 , 0.3 0 } are assigned with equal probability, that is, 0.5.For the case of binomial distribution, ten times Bernoulli trials with the probability 0.5 are adopted, and the integer values between 0 and 10 are linearly projected onto [−0.3 0 , 0.3 0 ].For the case of  distribution, the probability density function is defined as where (, ) is a  function, and the filter function  (0,1) () is defined as Three kinds of  distributions are considered, with  =  = 0.75,  =  = 1, and  =  = 4, respectively.The interval (0, 1) is linearly mapped to [−0.3 0 , 0.3 0 ].
Corresponding to these five distributions, the empirical cumulative distribution functions of relative errors are plotted in Figure 7. Comparing the ECDF curves at different effective correlation lengths, the same trend is observed as the Gaussian case.That is, the relative error decreases with the drop of effective correlation length.

Efficiency Analysis.
Using the same computing platform, the Monte Carlo method, the Neumann expansion method,  and the GNE method are employed to solve the same set of solutions (corresponding to 4340 samples), and the computation time is shown in Table 1.The number of random variables in the stochastic equation system differs depending on the effective correlation length.A fixed mesh is adopted in all of the computation such that the dimension of a single equation system is 434 × 434.
If the mesh density is doubled, then the dimension of each matrix is 1634 × 1634.The number of random variables remains the same in each case.The corresponding computation time is shown in Table 2.
The two data tables indicate that the computational efficiency of the GNE method (up to the third order) is much higher than Neumann expansion and Monte Carlo procedures.However, the memory requirement for higher order GNE schemes increases rapidly causing the computational efficiency to decrease.

Conclusions
An acceleration approach of Neumann expansion is proposed for solving stochastic linear equations.The new method,    The GNE method is equivalent to Neumann expansion for general stochastic equations, and for stochastic linear equations, both methods are also equivalent to the perturbation method.According to this equivalence relation, a rigorous error estimation is obtained, which is applicable to all the three methods.Numerical tests show that the third-order GNE expansion can provide good accuracy even for relatively large random fluctuations.

Figure 1 :
Figure 1: Cross section of the random material sample.

Figure 3 :Figure 4 :
Figure 3: The relation between the relative error and the spectral radius.

4 Figure 5 :
Figure 5: The relation between the relative error and the spectral radius.
and other B terms can be similarly formed.Matrices C  are determined by realignment of B  , B  , . . .and random variables   are formed by rearrangement of   1 ,   1  2 , . . .,   1  2 ⋅⋅⋅  .