Many problems in science and engineering involve spatially varying input data. Often, the input data is subject to uncertainties due to inherent randomness. For example, the details of spatial variations of properties and structure of engineering materials are typically obscure, and randomness and uncertainty are fundamental features of these complex physical systems. Under these circumstances, traditional deterministic models are rarely capable of properly handling this randomness and yielding accurate predictions. Therefore, in order to furnish accurate predictions, randomness must be incorporated directly into the model and the propagation of the resulting uncertainty between its input and output must be quantified accordingly. A model, particularly important in many applications, consists of the random input data in the form of a random field and a partial differential equation (PDE). The specific model problem considered here is
where is a nonlinear elliptic operator dependent on a random field
over a probability space, and is a solution. Examples of (1) include, but are not limited to, flow of water through random porous medium and modeling of the mechanical response of materials with random microstructure.
Over the past few decades, the critical need to seek solutions of (1
) has yielded a wealth of numerical approaches. Numerical methods for solving PDEs with random coefficients have been traditionally classified into three major categories: stochastic collocation (SC)babuvska2007stochastic ; nobile2008sparse ; nobile2008anisotropic , stochastic Galerkin (SG) ghanem2003stochastic ; babuska2004galerkin ; gunzburger2014stochastic ; xiu2002wiener ; xiu2003modeling and Monte Carlo (MC) babuska2004galerkin ; matthies2005galerkin ; kuo2016application . SC aims to first solve the deterministic counterpart of (1
) on a set of collocation points and then interpolate over the entire image space of the random element. Hence, the method is non-intrusive meaning that it can take advantage of existing legacy solvers developed for deterministic problems. Similarly, MC is non-intrusive as well since it relies on taking sample average over a set of deterministic solutions computed from a set of realizations of the random field. In contrast, SG is considered as intrusive since it requires construction of discretizations of both the stochastic space and physical space simultaneously and, as a result, it commonly tends to produce large systems of algebraic equations whose solutions are needed. However, these algebraic systems are considerably different from their deterministic counterparts and thus deterministic legacy solvers cannot be easily utilized.
We emphasize that all of the three categories discussed so far depend on the stochastic weak formulation of (1). In this work, alternatively, we take the variational viewpoint and reformulate problem (1) as a problem of seeking minimizers of the following functional
whose Euler-Lagrange equations coincide with (1) under suitable assumptions. Therefore, over an appropriate space, solving (1) is equivalent to minimizing (2). The use of variational formulations has become widespread in many areas of science and engineering due to their many advantages struwe1990variational ; reddy2017energy . First and foremost, the equations in the weak form are often applicable in situations when the strong form may no longer be valid. A case in point is modeling of microstructure evolution in materials, where fine scale oscillations may emerge, leading to highly irregular solutions ball1989fine ; muller1999variational . Second, variational formulations are known to be remarkably convenient for numerical computation as they often produce numerical methods capable of preserving, at least to some extent, the structure of the original problem simo2006computational ; marsden2001discrete .
In order to establish the existence of minimizers of in an appropriate space , by the direct method of calculus of variations dacorogna2007direct , it is sufficient to identify a minimizing sequence which satisfies two properties:
the sequence is compact under the weak topology on , i.e.,
This is often implied by the boundedness of the sequence (up to the extraction of a subsequence), i.e., for some constant independent of .
the functional is lower semicontinuous with respect to weak convergence, i.e.,
Given the above two properties, it is straightforward to verify that the function is indeed a minimizer of . The direct method is not only of theoretical importance. From the numerical point of view, it suggests that if we can identify a minimizing sequence each of which solves (2) over a finite dimensional subspace , i.e.,
then the above two properties ensure that converges weakly to the minimizer . That is, when is interpreted as a sequence of solutions to (2) over a sequence of finite dimensional spaces that approximates , the direct method of variational calculus automatically guarantees the numerical consistency. Bearing this in mind, numerical approximation to (2) boils down to minimizing over finite dimensional spaces with suitable optimization methods.
The fundamental difficulty in solving the stochastic optimization problem (2) is that the expectation often involves high dimensional integral which generally cannot be computed with high accuracy nemirovski2009robust . Thus, conventional nonlinear optimization techniques are seldom suitable for problems like (2
) since an inaccurate gradient estimation is usually detrimental to the convergence of the algorithms. In contrast, stochastic gradient descent (SGD) replaces the actual gradient by its noisy estimate, but is guaranteed to converge under mild conditionsbottou2010large ; bottou2018optimization ; kingma2014adam . The method can be traced back to the Robbins–Monro algorithm robbins1951stochastic
and has nowadays become one of the cornerstone for large-scale machine learningbottou2018optimization . However, due to the noisy nature of SGD iteration, a naive use of the algorithm in many instances suffers difficult tuning of parameters and extremely slow convergence rate nemirovski2009robust . In this article, we describe an application of SGD to construct numerical schemes for the solution of the variational stochastic problem (2). We also provide simple, yet powerful, strategies for efficient and robust SGD algorithms in the above context.
The reminder of the article is organized as follows. In Section 2, we setup the semilinear model problem and impose several running assumptions on the model. The variational reformulation of the model problem as a stochastic minimization problem is described in Section 3. Afterward, in Section 4, we propose to utilize the SGD to solve the minimization problem and discuss some useful technique for noise reduction and convergence acceleration for SGD. Finally, numerical benchmarks are presented in Section 5.
2 Model problem
We introduce a probability space where is the set of all events, is the -algebra consisting of all measurable events and a probability measure. We consider the following semilinear elliptic PDE with random coefficient defined in ,
where the domain is a bounded subset of , the boundary is either smooth or convex and piecewise smooth, the diffusion coefficient is a random field with continuous and bounded covariance functions and the nonlinear term is sufficiently smooth for almost surely all . We assume the solution to (3) exists and is unique.
We deal with the case of finite dimensional noise, i.e., there exists finitely many independent random variablessuch that
These random variables are often referred as the stochastic germ that bring randomness into the system. From the practical perspective, the finite dimensional noise assumption is reasonable since the random input often admits a parametrization in terms of finitely many random variables. From the theoretical perspective, by the Karhunen-Loève (KL) expansion mercer1909xvi : when the random field is square integrable with continuous covariance function, i.e.,
and the covariance function
is a well-defined continuous function of , then the random field can be approximated by the truncated KL expansion
where is the mean of the random field, are eigen-pairs of the covariance kernel
and are uncorrelated and identically distributed random variables with mean zero and unit variance.
by the Doob-Dynkin’s lemma oksendal2013stochastic . To define a suitable space of solution of the above problem, we introduce the physical space , i.e., the Sobolev space of functions with weak derivatives up to order and vanishing on the boundary. We also define the stochastic space , i.e., the space of -valued square integrable (with respect to ) random variables. The solution to (4
) is now thus defined in the tensor space
Now we make some technical assumptions on and . We denote the closure of .
is uniformly bounded and uniformly coercive, i.e., there exist constants such that
is uniformly bounded, i.e., there exists a constant such that
Furthermore, is uniformly Lipschitz continuous in , i.e., there exists a Lipschitz constant such that
Finally, is uniformly bounded from below, i.e., there exists a constant such that
As we shall see hereafter, these assumptions are crucial to guarantee the convergence of the SGD algorithm.
3 Direct method and polynomial chaos expansion
3.1 Direct method in calculus of variations
Following basic ideas of the calculus of variation, our starting point is to reformulate the stochastic PDE problem (4) as the minimization problem
with , and the expectation
taken with respect to the random vector. We first show that (5) has a minimizer and this minimizer satisfies the weak form of (4). The result is a simple application of the direct method in variational calculus dacorogna2007direct .
As we sketched in the introduction, it is sufficient to show that 1) there exists a minimizing sequence that converges weakly to some and 2) the functional is (weakly) lower semicontinuous dacorogna2007direct .
For weak convergence of the sequence , it suffices to uniformly bound under the norm. To this end, note that by Assumption 1 and the uniformly lower boundedness of
Integrating over and taking expectation of both sides lead to
Invoking the Poincaré’s inequality, there exist some constants and such that
Now note that is uniformly bounded (in ) since it is a minimizing sequence of (5), which implies is uniformly bounded and hence there exists such that
Next, we justify that is (weakly) lower semicontinuous. To this end, we define
and show both and are lower semicontinuous. To see the lower semicontinuity of , since converges weakly to in the Hilbert space , by definition (through choosing the test function )
Squaring both sides and then applying the Cauchy-Schwarz inequality lead to
If (the case when is trivial since ), the above inequality implies
i.e., is (weak) lower semicontinuous. Now for , by Taylor expansion in and the uniform boundedness of (Assumption 2), we have
The lower semicontinuity of follows immediately from the weak convergence of to . Therefore, is lower semicontinuous and hence has a minimizer .
It remains to be shown that satisfies the weak form (6). To this end, we consider the functional evaluated at for every , i.e., . A simple calculation shows that the Gateaux derivative satisfies
for some random variable . Since is uniformly bounded by Assumption 2, by the dominated convergence theorem
Owing to the fact that is a minimizer, the Gateaux derivative
which yields the weak form (6). ∎
Due to the infinite dimensionality of the solution space , it is not practical to solve (5) directly and hence we seek an approximate solution over a finite dimensional subspace , where is a generic index parameterizing the approximation accuracy. Consequently, the approximated minimization problem over the finite dimensional subspace becomes
Denote the minimizer of over the finite dimensional subspace . Since is a minimization sequence of (5), i.e.,
the sequence is compact on by Theorem 3.1, that is, converges weakly to , a minimizer of (5). That is to say, the consistency of the numerical approximation is guaranteed automatically under the framework of direct method of variational calculus. This serves as the theoretical foundation of the algorithm proposed in this work. Thus, the problem is reduced to finding a finite dimensional subspace minimizer as an approximation to the infinite dimensional space minimizer .
3.2 Polynomial chaos expansion
In this section, we make the above finite dimensional subspace approximation explicit. It is customary to construct approximations in the physical space by means of polynomials, in particular those with compact support, as is the case for the finite-element method. Undeniably, there exist various finite dimensional approximations to the stochastic space . For the sake of clarity, throughout this paper we adopt the generalized polynomial chaos (PC) expansion xiu2002wiener as a convenient approximation method in space . However, we emphasize that the general framework presented in this work extends naturally to other approaches for constructing approximants, such as piecewise polynomials expansions ghanem2003stochastic , multiwavelet decompositions le2004a ; le2004b to name a few.
The PC expansion is essentially a representation of second order random objects (e.g., random variables in or stochastic fields in ) xiu2002wiener ; xiu2003modeling ; le2010spectral . Given a random variable in , the PC expansion asserts that we can identify a set of -orthogonal univariate polynomial bases so that any function satisfying can be expressed as
in the sense, where the coefficients are
For instance, are Hermite polynomials when is normal and are Legendre polynomials when is uniform. The PC expansion can be generalized to the case of dimensional random vector with independent components. Specifically, for any function satisfying we can write
where are -variate polynomials involving products of those univariate polynomials associated with each component.
Now, given the dimensional random vector in (5), we can expand the random field as a generalized PC series
where the coefficient function
In practice, truncation of the PC series is required for numerical approximation. To this end, we define , the finite dimensional subspace spanned by the -orthogonal polynomials , i.e.,
Note that is determined by both the stochastic dimensionality and the highest order of the basis polynomials through
Hence, the dimensionality of the stochastic subspace can be very high when and are large. In order to further expand the coefficients , we approximate it over
the finite dimensional subspace of spanned by the bases . Therefore, the finite dimensional subspace of is over which we have a finite dimensional approximation
Note that in the notation of we omit the dependence on and in order to simplify the notation. For convenience, we define the vector valued function consisting of all bases of with the following numbering of index
Hence, the approximated solution can be written in the following compact form
where the coefficient (column) vector is
Over the finite dimensional subspace , the functional associated with becomes
Note that is indeed a function of , hence we rewrite as a function of the coefficients , denoted by . Therefore, minimizing is equivalent to minimizing the following function with respect to the coefficient vector ,
are the values associated with the linear part and the nonlinear part of (4), respectively. Here the gradient of with respect to is defined as
Finally, we make the following technical assumption regarding the basis functions.
For each , the physical-space basis satisfies
For each , the PC basis satisfies
The integrability conditions are satisfied in most settings. For example, when the physical domain is compact and are finite element bases, the integrability is readily verified. For the stochastic space, when is normal or uniform,
has finite moments of all orders.
4 Stochastic gradient descent for semilinear problem
4.1 Convergence of stochastic gradient descent
As mentioned in Section 3.1, in order to find an approximation to the solution of problem (3), it is sufficient to solve the stochastic optimization problem (9). The natural choice for optimizing the function is the stochastic gradient descent robbins1951stochastic , which is one of the most fundamental ingredients of large-scale machine learning bottou2010large ; bottou2018optimization . In this section, we discuss an application of SGD to the specific minimization problem (9
). For convenience, we denote the unbiased estimator ofby
so that . Instead of computing the deterministic gradient at each iteration, SGD simply requires the stochastic gradient for each iteration:
where is the learning rate. In the context of machine learning, each corresponds to a randomly picked example from the dataset. In the context of this article, we interpret as a realization of the stochastic germ . Note that the iterative sequence of coefficients is a sequence of random variables since each depends on for .
Intuitively, SGD works because, while each direction may not be one of the descent directions of , it is, however, a descent direction in expectation. It is clear that SGD is advantageous as it only requires computation of a single realization of the gradient at each iteration. Yet, it is a fundamental question whether SGD applied to the problem (9) produces a convergent sequence minimizing the function . To answer this question, we first present three important lemmas concerning the properties of the function and the gradient estimator .
Throughout the proof, denotes a generic positive constant that may differ by a scaling constant. Let be two arbitrary vectors of coefficients for the expansions over . For the linear part of , by Assumption 2,
where is the matrix -norm. To see that the matrix norm is finite, it is sufficient to bound terms of the form
Then, the function is strongly convex (in ), i.e., there exists a constant , such that for any
and hence has a unique minimizer .
Now, we are ready to present the main result concerning the convergence of SGD when applied for solving problem (9). Given Lemmas 4.1, 4.2 and 4.3, the proof of the following result is simply a straightforward application of Theorem in bottou2018optimization .
Some remarks about the above result are in order.
Diminishing learning rate has to be used to guarantee the convergence. The initial learning rate cannot be larger than a certain threshold.
We comment that the rate is the fastest convergence rate that the stochastic gradient descent can achieve agarwal2012information . However, the multiplicative constant can be improved by incorporating the second order information of the function . We will discuss more on this aspect in the next chapter.
A similar result can be obtained when the function is convex but not strongly convex. However, the convergence deteriorates to nemirovski2009robust .
It is well known that the SGD (12) suffers from the adverse effect of noisy gradient estimation. On the other hand, there is no particular reason to estimate the gradient only based on one realization of the random variables at each iteration. Therefore, it is natural to introduce a mini-batch at each iteration in order to “stabilize” the algorithm. That is, at the -th iteration we average the gradient over a batch of realizations of the random variable in order to obtain a less noisy gradient estimation
where and each is the -th realization in the mini-batch at the -th iteration of SGD. With the mini-batch gradient, we have the mini-batch SGD
Clearly, the mini-batch averaging reduces the variance of gradient estimation by a factor of but is also times more expensive than standard SGD. More sophisticated mini-batch strategies can be applied to accelerate SGD zhao2014accelerating ; dekel2012optimal . Our simulation results suggest that mini-batch is crucial for the convergence of SGD in order to take a relatively large learning rate.
4.2 The second order SGD
Further improvements to the SGD (12) may be achieved by way of incorporating the second order information pertaining to the function . Theorem 4.1 elucidates the fact that the constant appearing in the convergence rate depends on , which in turn depends on the condition number of the Hessian . This is similar to the deterministic optimization where the second order information is often incorporated to overcome the ill-conditioning of the optimization problem nocedal2006numerical . The same approach can be utilized in the stochastic setting by means of adaptively rescaling the stochastic gradients based on matrices capturing local curvature information of the function , so that the constant is significantly improved as a result. More precisely, we consider an iteration scheme
where is a symmetric positive definite approximation to the inverse of the Hessian . In fact, it was shown in bottou2005line that if is updated dynamically such that , then the multiplicative constant appearing in the convergence rate is independent of the condition number of the Hessian. It is true that approximation of is often based on a small set of samples and hence is very noisy. However, it has been long observed that the Hessian matrix need not be as accurate as the gradient in order to yield an effective iteration since the iteration is more tolerant to noise in the Hessian estimate than it is to noise in the gradient estimate. Therefore, it may be beneficial to incorporate partial Hessian information in the stochastic setting. To this end, we denote the unbiased estimator of the Hessian as
so that . Note that the estimator is indeed independent of whereas depends on in a nonlinear way. For convenience, we define two by matrices and with components
for . Then, we can express the Hessian estimator in the following block matrix form,