Numerical study in stochastic homogenization for elliptic PDEs: convergence rate in the size of representative volume elements

We describe the numerical scheme for the discretization and solution of 2D elliptic equations with strongly varying piecewise constant coefficients arising in the stochastic homogenization of multiscale composite materials. An efficient stiffness matrix generation scheme based on assembling the local Kronecker product matrices is introduced. The resulting large linear systems of equations are solved by the preconditioned CG iteration with a convergence rate that is independent of the grid size and the variation in jumping coefficients (contrast). Using this solver we numerically investigate the convergence of the Representative Volume Element (RVE) method in stochastic homogenization that extracts the effective behavior of the random coefficient field. Our numerical experiments confirm the asymptotic convergence rate of systematic error and standard deviation in the size of RVE rigorously established in [6]. The asymptotic behavior of covariances of the homogenized matrix in the form of a quartic tensor is also studied numerically. Our approach allows laptop computation of sufficiently large number of stochastic realizations even for large sizes of the RVE.


Introduction
Homogenization methods allow to derive the effective mechanical and physical properties of highly heterogeneous materials from the knowledge of the spatial distribution of their components [1,11]. In particular, stochastic homogenization via the Representative Volume Element (RVE) methods provide means for calculating the effective large-scale characteristics related to structural and geometric properties of random composites, by utilizing a possibly large number of probabilistic realizations [8,9,10,6,5]. The numerical investigation of the effective characteristics of random structures is a challenging problem since the underlying elliptic equation with randomly generated coefficients should be solved for many thousands realizations and for domains with substantial structural complexity to obtain sufficient statistics. Note that for every realization of the random medium one should construct a new stiffness matrix and right hand side to solve the discretized problem. Therefore, construction of fast solvers (which allow to confirm numerically the quantitative results in the stochastic homogenization) using the conventional computing facilities is a challenge.
This paper presents the numerical study of the stochastic homogenization of an elliptic system with randomly generated coefficients. Our approach is based on the FEM-Galerkin approximation of the 2D elliptic equations in a periodic setting by using fast assembling of the FEM stiffness matrix in a sparse matrix format, which is performed by agglomerating the Kronecker tensor products of simple 1D FEM discrete operators [12]. We use the product piecewise linear finite elements on the rectangular grid assuming that strongly varying piecewise constant equation coefficients are resolved on that grid. This scheme provides efficient approximation of equations with complicated jumping coefficients. The numerical analysis of the error in the Galerkin FEM approximation indicates the convergence rate O(h β ) in the L 2 -norm with 3/2 ≤ β ≤ 2.
The resulting large linear system of equations is solved by the preconditioned CG iteration, where the convergence rate is proven to be independent on the grid size and the relative variation in jumping coefficients, i.e., on the contrast. The preconditioned iterative solvers for the discrete elliptic systems of equations with variable coefficients have been long since in the focus in the literature on numerical methods for multi-dimensional and stochastic PDEs, see [13,14,16] and the literature therein.
In this paper, we consider an ensemble of two-valued random coefficient fields, which is based on independently and uniformly placed (and thus overlapping) axi-parallel square inclusions of fixed side length. We investigate the RVE method that (approximately) extracts the effective (i.e. large-scale) behavior of the medium in form of the deterministic and homogeneous matrix A hom from a given (stationary and ergodic) ensemble. This method produces an approximation to A hom by solving D = 2 elliptic equations on a torus of (lateral) size L with a specific right hand side (the corrector equation), by taking the spatial average of the flux of these solutions, and by taking the empirical mean over N independent realizations of this coefficient field under the naturally periodized version of the ensemble. This is an approximation in so far as the outcome is still random (as quantified by the standard deviation of the outcome of a single realization) and that the periodic boundary conditions affect the statistics (which we call the systematic error). In [6], Gloria, Neukamm and the last author rigorously derived upper bounds how the standard deviation and the systematic error decrease with increasing RVE size L. Our numerical experiments confirm the scaling of these bounds. Since numerically, there is no access to exact values of the variance (or standard expectation) or the expectation, we replace these quantities by their empirical counterparts for a large number of realizations N . We thus first provide numerical evidence that these quantities have saturated in N , and second that their limiting values display the predicted scaling in L.
In work [4] by Duerinckx, Gloria and the last author, it was worked out that the properly rescaled variance of the output of the RVE converges as L ↑ ∞ to a quartic tensor Q that governs the leading-order fluctuations of any solution. In this paper, we show how the symmetry properties of the ensemble yields symmetry properties of Q (and its approximation). Also, a convergence rate was rigorously established in that work, and is being numerically investigated here. While the number of samples N considered here seems too low to reach saturation in the empirical variance, the findings are at least not in contradiction with the theoretic results.
In numerical tests on the stochastic properties of the 2D RVE method we study the asymptotics of empirical variance versus the size of RVE L ≤ 128, and of the systematic error versus the number of realizations N up to N = 10 5 . Furthermore, we estimate the convergence of the quartic tensor by implementing a large number of stochastic realizations. The proposed techniques allow to compute a sufficiently large number of realizations of random coefficient fields with a large number of overlapping inclusions up to L 2 = 128 2 corresponding to the stiffness matrix size 513 2 × 513 2 using MATLAB on a moderate computer cluster.
The numerical investigation of the stochastic homogenization problem attracts interest and becomes an active field of research, see the survey [1] and references therein. Recently the numerical solution of the corrector-type problem, in the context of homogenization of the diffusion equation with spherical inclusions by using boundary element methods and the fast multipole techniques has been considered in [3].
The rest of the paper is organized as follows. In Section 2, we address the problem setting and define the elliptic equations of stochastic homogenization. Section 3 describes the Galerkin-FEM discretization scheme based on the fast matrix generation by using sums of Kronecker products of single-dimensional matrices. We also outline the PCG iteration applied in the computer simulations and provide numerics on the FEM discretization error. Section 4 introduces the computational scheme for the stochastic average coefficient matrix. Furthermore, in Section 4.3 we describe the construction and properties of the covariances of the homogenized matrix in the form of a quartics tensor. Section 5 presents results of numerical experiments on the empirical average and systematic error at the limit of a large number of stochastic realizations. The asymptotic of the quartic tensor versus the leading order variances is analyzed numerically in Section 5.2. Conclusions outline the main results of the paper. In Appendix the spectral properties of the discrete stochastic elliptic operators are studied numerically for a sequence of stochastic realizations.

Elliptic equations in stochastic homogenization
In this section, we describe the problem setting in the stochastic homogenization theory. For given f ∈ L 2 (Ω) such that Ω f (x)dx = 0, we consider the class of model elliptic boundary value problems on the d-dimensional torus Ω : In this paper we focus on the special class of elliptic problems (2.1) arising in stochastic homogenization theory where the highly varying coefficient matrix and the right-hand side are defined by a sequence of stochastic realizations as described in [8,9,10,6,5], see details in Sections 3 and 4.
In what follows, we present the numerical analysis for 2D stochastic homogenization problems (2.1) with periodic boundary conditions on Γ = ∂Ω, in the form where the scalar function a(x) > 0 is piecewise constant in Ω. The efficient numerical simulation pre-supposes the fast numerical solution of the equation (2.1) in the case of many different coefficients A(x) and right-hand sides, generated by certain stochastic procedure, and in the calculation of various functionals on the sequence of solutions φ. In this problem setting the bottleneck task is fast and accurate generation of the FEM stiffness matrix in the sparse matrix form, which should be re-calculated many hundred if not thousand times in the course of stochastic realizations.
In asymptotic analysis of stochastic homogenization problems the coefficient and the right-hand side are chosen in a specific way, see [10,8] for the particular problem setting. In this paper, we describe the conventional 2D FEM discretization scheme in the rescaled domain Ω = [0, 1) 2 . Given the size of representative volume elements L = 1, 2, 3, . . . we generate (possibly overlapping) decomposition of the domain Ω into L 2 equal unit cells G s , s = 1, . . . , L 2 , each of size 2α L × 2α L , whose centers are distributed randomly (by the Poisson distribution) within the supercell Ω, taking into account periodicity in both spacial variables x 1 and x 2 . Here the overlap factor satisfies α ≤ 1/2. Stochastic characteristics of the system can be estimated at the limit of a large number L.
We consider a sequence of random coefficient distributions {G s } n , numbered by n = 1, . . . , N , where the particular set {G s } = {G s } n for fixed n will be called a realization. For any fixed realization define the covered domain and the respective coefficient 3) The stochastic model is specified by the choice of the overlap constant α ∈ (0, 1/2] and the scaling factor λ ∈ (0, 1]. In the following, the constant λ will be fixed in the interval 0.1 ≤ λ ≤ 0.8. We denote the "stochastic" elliptic operator for the particular realization by where the corresponding 2 × 2 coefficient matrix A (n) (x) is defined by (2.4) and the diagonal matrix coefficient takes the form (2.5) We use the notation A (n) (x) for the "stochastic part" of a matrix associated with the diagonal coefficient a (n) (x), i.e.
Now the elliptic equations in stochastic homogenization are formulated as follows. Fixed coefficient A (n) (x), for i = 1, 2 solve the periodic elliptic problems in Ω, where the directional unit vectors e i , i = 1, 2, are given by e 1 = (1, 0) T and e 2 = (0, 1) T , see Section 4 for more details. The variational form of the equation (2.6) reads as follows The equation (2.6) can be also written in the classical form (2.1) The problem setting remains verbatim in the d-dimensional case, d > 2. In this case the equation (2.8) takes the same form, where a d × d coefficient matrix is given by and e i ∈ R d , i = 1, . . . , d, represent the set of directional unit vectors in R d .

Matrix generation and iterative solution
In this section we describe the FEM discretization scheme and the fast matrix generation approach based on the use of tensor Kronecker products of "univariate" matrices.

Galerkin FEM discretization
First, we introduce the uniform m s × m s rectangular grid Ω hs in Ω with the grid size h s = 1 ms−1 , such that m s = m 0 L + 1, m 0 = 2 p 0 , i.e. h s = 1 m 0 L . We assume that the unit cell G s , s = 1, . . . L 2 , of size 2α L × 2α L adjust the square grid Ω hs , such that the center c s In the following numerical examples we normally use the overlap constant α = 1/4. For α = 1/2, the maximal size of the unit cell is given by 1 L × 1 L , which contains m 0 + 1 grid points in each spacial direction leading to m 1 × m 1 rectangular grid with m 1 = m 0 L + 1.
The FEM discretization of the elliptic PDE in (2.8) can be constructed, in general, on the finer grid Ω h compared with Ω s , which serves for the resolution of jumping coefficients. To that end, we introduce the m 1 × m 1 rectangular grid Ω h with the mesh size h = 1 m 1 −1 , m 1 ≥ m s , that is obtained by a dyadic refinement of the grid Ω s , such that the relation holds, implying h s = 2 p h. Now the grid-size of the unit cell G s on the finer grid Ω h is given by (m 0 2 p + 1) × (m 0 2 p + 1). Given a finite dimensional space X ⊂ H 1 (Ω) of tensor product piecewise linear finite elements X = span{ψ µ (x)} associated with the grid Ω h , with µ = 1, ..., M d , M d = m d 1 , for d = 2 incorporating periodic boundary conditions, we are looking for the traditional FEM Galerkin approximation of the exact solution in the form where u = (u 1 , . . . , u M d ) T ∈ R M d is the unknown coefficients vector. Fixed realization of the coefficient a (n) (x), for i = 1, 2 we define the Galerkin-FEM discretization with respect to X of the variational equation (2.7) by where the Galerkin-FEM matrix A generated by the equation coefficient A (n) (x) is calculated by using the associated bilinear form and Corresponding to (2.8) and (3.3), we represent the stiffness matrix A in the additive form where A ∆ represents the M d × M d FEM Laplacian matrix in periodic setting that has the standard two-terms Kronecker product form. Here matrix A s provides the FEM approximation to the "stochastic part" in the elliptic operator corresponding to the coefficient a (n) (x), see (2.4). The latter is determined by the sequence of random coefficients distribution in the course of stochastic realizations, numbered by n = 1, . . . , N .
In the case of complicated jumping coefficients the stiffness matrix generation in the elliptic FEM usually constitutes the dominating part of the overall solution cost. In the course of stochastic realizations the equation (3.2) is to be solved many hundred or even thousand times, so that every time one has to update the stiffness matrix A and the righthand side f .
Our discretization scheme computes all matrix entries at the low cost by assembling of the local Kronecker product matrices obtained by representation of a (n) (x) as a sum of separable functions. This allows to store the resultant stiffness matrix in the sparse matrix format. Such a construction only includes the pre-computing of tri-diagonal matrices representing 1D elliptic operators with jumping coefficients in periodic setting. In the next sections, we shall describe the efficient construction of the "stochastic" term A s .

Matrix generation by using Kronecker product sums
To enhance the time consuming matrix assembling process we apply the FEM Galerkin discretization (3.3) of equation (2.8) by means of the tensor-product piecewise linear finite elements where ψ µ (x ) are the univariate piecewise linear hat functions 1 . The M d × M d stiffness matrix is constructed by the standard mapping of the multi-index µ into the long univariate index 1 ≤ µ ≤ M d for the active degrees of freedom in periodic setting. For instance, we use the so-called big-endian convention for d = 3 and d = 2 respectively. In what follows, we consider the case d = 2 in more detail.
In our discretization scheme we calculate the stiffness matrix by assembling of the local Kronecker product terms by using representation of the "stochastic part" in the coefficient a (n) (x) as an R-term sum of separable functions. To that end, let us assume for the moment that the scalar diffusion coefficient a(x 1 , x 2 ) can be represented in the separate form (rank-1 representation) Then the entries of the Galerkin stiffness matrix A = [a µν ] ∈ R M d ×M d can be represented by which leads to the rank-2 Kronecker product representation where ⊗ denotes the conventional Kronecker product of matrices, see Definition 3.1 below.
Here A 1 = [a µ 1 ν 1 ] ∈ R m 1 ×m 1 and A 2 = [a µ 2 ν 2 ] ∈ R m 2 ×m 2 denote the univariate stiffness matrices and S 1 = [s µ 1 ν 1 ] ∈ R m 1 ×m 1 and S 2 = [s µ 2 ν 2 ] ∈ R m 2 ×m 2 define the weighted mass matrices, for example Definition 3.1 Recall that given p 1 × q 1 matrix A and p 2 × q 2 matrix B, their Kronecker product is defined as a p 1 p 2 × q 1 q 2 matrix C via the block representation Let us discuss in more detail the calculation of the 1D stiffness matrices A 1 and A 2 in the case of variable 1D coefficients. We choose the Galerkin FEM with m = m 1 piecewise-linear hat functions {ψ µ 1 } in periodic setting in Ω = [0, 1), constructed on a uniform grid with a step size h = 1/m, and nodes x µ 1 = h µ 1 , µ 1 = 1, . . . , m. If we denote the diffusion coefficient by a(x 1 ), then the entries of the exact stiffness matrix A 1 read as We assume that the coefficient remains constant at each spatial interval [x µ 1 −1 , x µ 1 ], which corresponds to the evaluation of the scalar product above via the midpoint quadrature rule yielding the approximation order O(h 2 ). Introducing the coefficient vector a = [a µ 1 ] ∈ R m , a µ 1 = a(x µ 1 −1/2 ), µ 1 = 1, . . . , m, where x i 1 −1/2 is the middle point of the integration interval, the symmetric tridiagonal matrix of interest can be represented by By simple algebraic transformations (e.g. by lamping of the mass matrices) the matrix A can be simplified to the form (without loss of approximation order) where D 1 , D 2 are the diagonal matrices with positive entries. This representation applies in particular to the periodic Laplacian.
In the general case, the piecewise constant stochastic coefficient can be represented as an R-term sum of separable coefficients. This leads to the linear system of equations The representation in (3.7) can be further simplified to the anisotropic Laplacian type matrix which will be used as a prototype preconditioner for solving the target linear system (3.8).
Taking into account the rectangular structure of the grid, we use the simple finitedifference (FD) scheme for the matrix representation of the Laplacian operator −∆. In this case the scaled discrete Laplacian incorporating periodic boundary conditions takes the form such that the entries of the "periodization" matrix P (1) ∈ R m 1 ×m 1 are all zeros except Here I m 1 ∈ R m 1 ×m 1 is the identity matrix, ∆ 1 = ∆ 2 is the 1D finite difference Laplacian (endorsed with the Neumann boundary conditions), and ⊗ denotes the Kronecker product of matrices, see Definition 3.1. We say that the Kronecker rank of both A in (3.7) and A ∆ in (3.9) equals to 2.
Notice that the m 1 × m 1 Laplacian matrices for the Neumann and periodic boundary conditions in 1D read as respectively.
In the d-dimensional case we have the similar Kronecker rank-d representations. For example, in the case d = 3 the "periodic" such that its Kronecker rank equals to 3, and similar for the arbitrary d ≥ 3.

Fast matrix assembling for the stochastic part
The Kronecker form representation of the "stochastic" term in (3.3) further denoted by A s is more involved. For given stochastically chosen distribution of overlapping cells G s , s = 1, . . . , L 2 , we construct the minimal non-overlapping decomposition of the full covered grid domain G = ∪ L 2 s=1 G s colored by gray in Figure 2.
Here m 0 = 2 p + 1, and p = 0, 1, 2, . . ., is fixed as above by relation m 1 − 1 = (m s − 1)2 p . In this construction, the non-overlapping elementary cells S k for different k are allowed to have the only common edges of size m 0 . Notice that in the case of non-overlapping decomposition (2.2) the set of cells {S k } coincides with the initial set {G s } which allows to maximize the size m 0 × m 0 of each S k , k = 1, . . . , L 2 , to the largest possible, i.e. to m 0 = m 0 2 p + 1.
To finalize the matrix generation procedure for A s , we define the local m 0 × m 0 matrices representing the discrete Laplacian with Neumann boundary conditions,

and the diagonal matrix
Let the subdomain S k be supported by the index set I of size m 0 × m 0 for k = 1, . . . , K. Introduce the block-diagonal matrices Q k ∈ R m 1 ×m 1 and I k ∈ R m 1 ×m 1 by inserting matrices Q m 0 and I m 0 as diagonal blocks into m 1 × m 1 zero matrix in the positions I (1) k , respectively. Now the stiffness matrix A s is represented in the form of a Kronecker product sum as follows, where is the "periodization" matrix in 2D. In a d-dimensional case the representation (3.13) generalizes to a sum of d-factor Kronecker products where P (d) is the "periodization" matrix in d dimensions, constructed as the d-term Kronecker sum similar to the case d = 2.
The Kronecker product form of (3.9) and (3.13) leads to the corresponding Kronecker sum representation for the total stiffness matrix A. This allows an efficient implementation of the matrix assembly and low storage for the stiffness matrix preserving the Kronecker sparsity. Hence, it proves the following storage complexity for the matrix A. Here, in general, the number K of elementary cells 2 is larger than L 2 , and it coincides with L 2 only in the case of non-overlapping decomposition G = ∪ L 2 s=1 G s , where different patches G s are allowed to have joint pieces of boundary but no overlapping area.
In the general case d ≥ 2 and K ≥ L d , the Kronecker rank of the matrix A is bounded by rank Kron (A) ≤ d K.

Numerical analysis of the FEM approximation error
We tested convergence of the solutions on a sequence of dyadic refined grids, for the fixed configuration of coefficients and the right hand side given by f = sin (2πx)

Preconditioned CG iteration
Let the right-hand side in (3.3) satisfy F, 1 = 0, then for a fixed m, the equation has the unique solution. We solve this equation by the preconditioned conjugate gradient (PCG) iteration (routine pcg in Matlab library) with the preconditioner where δ > 0 is a small regularization parameter introduced only for stability reasons (can be ignored in the theory) and I is the M d × M d identity matrix. It can be proven that the condition number of preconditioned matrix is uniformly bounded in m 1 , L and in the number of stochastic realizations n = 1, . . . , N . The particular estimates on the condition number in terms of a parameter λ can be derived by introducing the average coefficient where a + (x) and a − (x) are chosen as majorants and minorants of a (n) (x) in (2.4), respectively. The following simple result holds.
Proof. Lemma 4.1 in [15] shows that the preconditioner A 0 generated by the coefficient a 0 (x) = 1 2 (a + (x) + a − (x)) allows the condition number estimate The preconditioner B corresponds to the choice a + (x) = 1 and a − (x) = λ, hence, we obtain a 0 (x) = 1+λ 2 and the result follows. The PCG solver for the system of equations (3.8) with the shifted discrete Laplacian as the preconditioner demonstrates robust convergence with the rate q 1. In the practically interesting case α ≈ 0.5 we found that q does not depend on λ. This can be explained by the fact that in this case the total overlap in all subdomains covers the large portion of the computational box Ω. In all numerical examples considered so far the number of PCG iterations was smaller than 10 for the residual stopping criteria δ = 10 −8 . We use the univariate grid size m 1 = m s , corresponding to the choice p = 0 in (3.1) which is fine enough to resolve geometry for larger L.

Asymptotic convergence to the stochastic average
In this section, we describe the computational scheme for calculation of the homogenized coefficient matrix for each stochastic realization.

Computational scheme for the stochastic average
For fixed stochastic realizations specifying the variable part in the 2 × 2 coefficient matrix A (n) (x), n = 1, . . . , N , we consider the problems Taking into account (2.4), where the diagonal of A (n) (x) is defined in terms of the scalar function a (n) (x), we arrive at   ij ] ∈ R 2×2 , i, j = 1, 2, with the constant entries is given bȳ which implies the representation for matrix elements The latter leads to the entry-wise representation of the matrixĀ (n) = [ā The representation (4.3) ensures the symmetry of the homogenized matrixĀ (n) , i.e.ā ji . Indeed, we calculate the difference between the scalar product of the first equation in (4.1) with φ 2 , and the second equation in (4.1) with φ 1 , which then implies the desired property via integration by parts, and taking into account the relation (2.5), In numerical implementation, we apply the Galerkin scheme for FEM discretization of equation (4.1) its right-hand side. We use the same quadrature rule for computation of integrals in (4.3) thus preserving the symmetry in the matrix A (n) inherited from the exact variational formulation (see argument above and Section 4.3 for the more detailed discussion).
Integrals over Ω in (4.2), (4.3) for the matrix entries (Ā (n) ) i,j , i, j = 1, 2, are calculated (approximately) by the scalar product of the N -vector of all-ones with the discrete representation of integrand on the grid Ω h , see  To complete this section, we check numerically that the FEM discretization scheme preserves the symmetry in the matrixĀ (n) L for fixed L if the discrete system of equations (3.2) is solved accurately enough. Table 4.1 demonstrates that the symmetry in the matrixĀ (n) L with fixed n is recovered on the level of residual stopping criteria δ > 0 in the preconditioned iteration for solving the discrete system of equations. For this calculation we set L = 4, m 0 = 8, α = 0.5 and λ = 0.2.

Asymptotic of systematic error and standard deviation
The set of numerical approximations {Ā  In what follows, we use the entry-wise notation for d × d matrices A = [a ij ], i, j = 1, . . . , d, for example, Ā L = [ā L,ij ] andĀ L,ij ], etc. In terms of square expectations, the convergence rate for the computable quantities can be estimated by, see [6], We numerically study the asymptotic of both terms on the right-hand side of (4.7) separately by considering the random part of the error, 8) and the systematic error where Ā L L is calculated for large enough N .

Covariances of the homogenized matrix in the form of quartic tensor
Let · L be an ensemble of uniformly elliptic symmetric coefficient fields on the d-dimensional torus [0, L) d . Assume that it is invariant under translation (stationary) and under the group G of all orthogonal transformations R of R d that leave the (hyper-)cube [0, L) d invariant (this is generated by rotations in one of the Cartesian two-dimensional planes and reflections along any Cartesian hyper-plane) in the sense of (4.22) below. In case of isotropic (i.e. scalar) coefficient fields A(x), (4.22) turns into A(R·) and A have the same distribution under · L , which is certainly the case for the ensembles we consider numerically. Let X be a finite-dimensional space of functions on the torus [0, L) d of side-length L with square-integrable gradients, e.g. coming from continuous, piecewise affine Finite Elements. For a given realization A(x) = A (n) (x) of the coefficient field and any direction i = 1, · · · , d, where e i denotes the unit vector in direction i. If X contains the constant functions (as would be the case for the Finite Element space), φ i has to be normalized to be unique, e.g. Here it is important to treat periodicity correctly: In practice, one identifies functions on [0, L) d with functions on R d that are periodic in each (Cartesian) argument of period L, hence if the node α is such that one of the adjacent triangles crosses the boundary of the periodic cell [0, L) d ⊂ R d , then there is a piece of φ α that appears on the other side. If a quadrature rule is used for computing the stiffness matrix, it is important that the same one is used for approximation of the right-hand side. Let us consider the d × d matrixĀ L = [ā L,ij ] =Ā L (A) defined through (see also (4.2)) (where again, the same quadrature rule should be used). Then we have for every realization Let us consider the ensemble average Ā L L , which by the law of large numbers is given by (see also (4.5)) (4.14) almost surely, whereĀ (n) L come via (4.12) from independent realizations A = A (n) according to the distribution · L . Suppose that the finite-dimensional space X is invariant under reflections in the coordinate directions in the sense of (4.23) below. This imposes a more serious restriction on the Finite Element space, namely that it is based on a subdivision of the torus [0, L) d into axi-parallel cubes (instead of triangles) and that the function space on each cube is spanned by functions that are multi-linear in the Cartesian coordinates (as opposed to affine). If this condition is satisfied, then we have Ā L L is isotropic, i.e. ā L,ij L = λ L δ ij (4.15) for some λ L ∈ (0, L).
We are interested in the covariances of the entries ofĀ L , and note that by the law of large numbers More precisely, we are interested in its rescaled version which is easier to understand as the four-linear form We claim that it has the invariance propertȳ Q L (Rη, Rξ, Rη , Rξ ) =Q L (η, ξ, η , ξ ).
Argument for (4.15). Identifying the points on the torus with [0, L) d ⊂ R d , let G denote the subgroup of the orthogonal group that leaves [0, L) d invariant. According to our assumption, for any R ∈ G, R t A(R·)R and A have the same distribution under · L , According to our assumption on X we have For a fixed vector ξ ∈ R d , we consider φ ξ := ξ i φ i (Einstein's summation convention) and note that in view of (4.10), for given realization A = A (n) , the function φ ξ = φ ξ (A) (at least up to additive constants) is characterized by We now argue that φ transforms under R ∈ G as follows Indeed, this relies on the straightforward orthogonal transformation rule According to (4.23) and (4.24) (with ξ replaced by Rξ) the left-hand side vanishes for all ψ ∈ X; hence by the characterization (4.24) applied to the right-hand side, we obtain (4.25). We now argue note that from (4.25) we obtain for the gradient ∇φ Rξ (A; Rx) = ∇φ ξ (R t A(R·)R; x) and thus for the flux q ξ (A; x) := A(ξ + ∇φ ξ (A; x)) the transformation rule from which we obtain by definition (4.12) that (4.26) According to (4.22) this yields the following invariance property for the symmetric matrix Ā L L Ā L L Rξ = R Ā L L ξ.
Since this holds for all ξ ∈ R d and all R ∈ G, by an argument of elementary algebra, we obtain the isotropy of Ā L L , cf (4.15).
Argument for (4.16). This follows from (4.26)  The identity (4.20) follows from the symmetry of the covariance in its two arguments. The vanishing of the eight entries stated in (4.17) and (4.18) follows from (4.16) applied to the reflection R ∈ G given by Re 1 = −e 1 and Re 2 = e 2 . The identity (4.21) follows from (4.16) applied to the reflection R ∈ G given by Re 1 = e 2 and Re 2 = e 1 .

Numerical study of stochastic homogenization
In this section, we estimate numerically the mean constant coefficient in the system (2.8) depending on L and other model parameters at the limit of N → ∞, see [10,8] for the respective problem setting.
Recall that the homogenization problem is solved in the unit square Ω =     Table 5.1: CPU times (sec) versus the number of inclusions (i.e., L 2 ) for generating the stiffness matrix, the right-hand side, and for the solution of the discretized system for the case of overlapping inclusions. Tolerance ε = 10 −8 . Table 5.1 presents the CPU times for generating the stiffness matrix, the right-hand side (RHS), and for the solution of the discretized system for the case of overlapping inclusions, for tolerance ε = 10 −8 . Number of inclusions (L 2 ) varies from 16 to 16384. The latter is computed on a mesh of size 513 × 513. We observe that matrix generation takes the dominating time.

Systematic error and empirical variance versus L
In what follows, we numerically check the theoretical convergence rate (4.7), in form of checking (4.8) and (4.9) separately. Figure 5.2 serves to illustrate the asymptotic convergence of the systematic error see (4.9), at the limit of large L. Since we do not have access to the ensemble averages Ā L L , we take empirical averagesĀ N L for large enough N , (cf. (4.6)) as a proxy. Furthermore, due to the fact that A hom is not computable we compare the differences in ā L,11 L computed on a sequence of increasing values of L.

The asymptotic of quartic tensor vs. leading order variances
In this section, we consider the convergence of the quartic tensorQ L , representing covariances of the matrixĀ L , to its leading order variances Q hom , see Section 4.3. For the large number of realizations N , the computable approximation,Q N L ∈ R 2×2×2×2 , to the scaled quartic tensor is defined bȳ so that by the central limit theoremQ The equivalent matrix representation ofQ N L is obtained by setting the operation ⊗2 in (5.3) as the Kronecker product of matrices (see Definition 3.1), further denoted bȳ . . , 4. In our numerical tests we shall check the asymptotic behavior which can be expected at the limit of large size L of the RVE, see [6,4]. It is worth to note that the quartic tensorQ L can be calculated at no further cost than the effective homogenized matrixĀ L .

Conclusions
We present the numerical scheme for discretization and solution of 2D elliptic equations with strongly varying piecewise constant coefficients arising in stochastic homogenization of multiscale composite materials. The resulting large linear system of equations is solved by the preconditioned CG iteration with the convergence rate that is independent of the grid size and of the variation in jumping coefficients. For a fixed size of the representative volume element, our approach allows to avoid the generation of the new FEM space in each stochastic realization. For every realization, fast assembling of the FEM stiffness matrix is performed by agglomerating the Kronecker tensor products of 1D FEM discretization matrices. The resultant stiffness matrix is maintained in a sparse matrix format.
Our numerical scheme allows to investigate the asymptotic convergence rate of significant quantities of stochastic homogenization process in the course of a large number of realizations (of the order of N = 10 5 ) and for large sizes of the representative volume elements up to L = 128, corresponding to the number of inclusions 16384 and matrix size 513 2 × 513 2 . Note that for every realization a new matrix generation and solution of the respective linear system is performed.
Our numerical experiments study the asymptotic convergence rate of systematic error and standard deviation in the size of RVE, rigorously established in [6]. In particular, we confirm in various numerical tests the theoretical asymptotic estimates, see Section 4.2, concerning the convergence rate O(1/L) for the empirical variance at the limit of large L, but with a moderate number of stochastic realizations N , and the asymptotic CL −2 ln 2 L in the case of large N .
The asymptotic behavior of covariances of the homogenized matrix in the form of quartic tensor are studied numerically. In particular, we consider the asymptotic of the quartic tensor versus the leading order variances, computed for the large number of stochastic realizations up to N = 10 4 . In this way, the asymptotic O(L −d ln d L), for d = 2, is confirmed on a sequence of increasing sizes of the RVE, up to L = 64.
The stochastic characteristics of the system are analyzed for a range of intrinsic model parameters like the number of realizations, the size of periodic representative volume element, the jump-ratio in the stochastic equation coefficients (contrast) and various grid discretization parameters. The presented numerical scheme allows to perform large scale simulations using MATLAB on a moderate computer cluster. The tensor-based numerical techniques to matrix generation presented in this paper can be extended to 3D and higher dimensional problems.
7 Appendix: spectral density of a stochastic operator Spectral properties of the randomly generated family of elliptic operators are important in many applications, in particular, in stochastic homogenization of time dependent PDEs.
In what follows, we analyze numerically the average behavior of the so-called density of spectrum (DOS) for the family of stochastically generated 2D elliptic operators {A n } for the large sequence of stochastic realizations n = 1, . . . , N . The DOS provides the important spectral characteristics of the stochastic differential operator which accumulate the significant information on the static and dynamical characteristics of the complex physical or molecular system. Here we numerically demonstrate the convergence of DOS to the sample average function at the limit of large number of stochastic realizations with fixed L. Our second goal is the numerical study of the DOS depending on the increasing number L. We use the simple definition of DOS for symmetric matrices, see [18,2], where δ is the Dirac function and the λ j 's are the eigenvalues of A = A T ∈ R m×m , Au j = λ j u j , j = 1, . . . , m, assumed to be labeled non-decreasingly.
In the presented analysis we employ the commonly used class [18] of blurring approximations to the spectral density φ(t) by using regularization via a Gaussian function with width parameter η > 0, where the choice of small regularization parameter η depends on the particular problem setting. The DOS in (7.1) will be approximated by Gaussian broadening, φ(t) → φ η (t) := 1 m  Finally, we notice that the specific feature of the homogenized DOS is that the sample average differs substantially from DOS for the operator with the homogenized coefficient.