Calculating the Malliavin derivative of some stochastic mechanics problems

The Malliavin calculus is an extension of the classical calculus of variations from deterministic functions to stochastic processes. In this paper we aim to show in a practical and didactic way how to calculate the Malliavin derivative, the derivative of the expectation of a quantity of interest of a model with respect to its underlying stochastic parameters, for four problems found in mechanics. The non-intrusive approach uses the Malliavin Weight Sampling (MWS) method in conjunction with a standard Monte Carlo method. The models are expressed as ODEs or PDEs and discretised using the finite difference or finite element methods. Specifically, we consider stochastic extensions of; a 1D Kelvin-Voigt viscoelastic model discretised with finite differences, a 1D linear elastic bar, a hyperelastic bar undergoing buckling, and incompressible Navier-Stokes flow around a cylinder, all discretised with finite elements. A further contribution of this paper is an extension of the MWS method to the more difficult case of non-Gaussian random variables and the calculation of second-order derivatives. We provide open-source code for the numerical examples in this paper.


Introduction
The classical derivative is a fundamental tool of calculus that is widely used across every field of mathematics, science and engineering. Various generalisations and extensions of the classical derivative, e.g. local and/or partial Frechét and Gâteaux derivatives [1], are now common tools in the repertoire of practitioners working in many fields. Modern extensions such as fractional and non-local derivatives are finding increasing use in several fields of science and technology, see e.g. [2][3][4][5][6]. The semi-inverse method of [7] is a powerful tool for the establishment of variational principles (Euler-Lagrange) from governing equations for physical problems.
By contrast, the Malliavin calculus [8], an extension of the notions of classical calculus of variations to stochastic processes, is certainly less widely known. In our view, this is probably because the vast majority of papers written on the subject require study of mathematics and stochastics to an advanced level. However, we think that Malliavin calculus deserves a wider audience. The objective of this paper then is introduce the Malliavin derivative as a useful numerical tool for practitioners to understand the behaviour of stochastic PDEs in mechanics, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 In a practical way, if m is a random variable with probability density function f m , Eq (1) writes: As we will see, the Malliavin Weight Sampling method (MWS) [16] allows the evaluation of the sensitivity of the expected value of the quantity of interest with respect to the mean value of the stochastic parameter m as [9,11,16,18]: where q m is the Malliavin weight for the parameter m and " m is the mean of m. Under certain condition of regularity [11,20,21] when the probability density function (PDF) of the parameter m is known, the Malliavin weight q m associated can be computed directly from the PDF of m. This approach can be viewed as an integration by parts, and is a direct result of Malliavin calculus where we take the derivative of random functions rather than the classical derivative. We emphasise again the quite different nature of the above derivative Eq (3) to the classical notion of a derivative from elementary calculus.
In Eq (3), we suppose that the quantity of interest J does not depend explicitly on the parameter m. Later we introduce a more general equation Eq (35) that must be considered if in fact J does depend on m.
The simplest approach to calculate E½Jq m , and the one we use exclusively in this paper, is to use the standard Monte Carlo estimator; that is, take Z independent and identically distributed (iid) realisations m z of m, solve for J z ≔ J(u(m z )) before taking the sample mean of the set of realisations {J 1 , . . ., J Z }: Jðm z Þ Á q m ðm z Þ: ð4Þ From the central limit theorem, the error in Eq (4) is normally distributed with variance Z −1 V where V is the variance of Jq m .
What will not be clear to the reader at this stage is how to determine the Malliavin weights. Through a simple practical examples in the next section, we will explain how to use the MWS method, determine the specification of the weights for both Gaussian and non-Gaussian distributions on the parameter m, and thus calculate the Malliavin derivative Eq (3).

Kelvin-Voigt model
The Kelvin-Voigt constitutive model with Young's modulus E, viscosity η and loading stress σ can be written as the following linear ordinary differential equation: A schematic of this model is shown in Fig 1. The initial condition on the strain is (t = 0) = 0 and we study the response of the system for time t 2 [0, T]. Our quantity of interest functional is the value of the strain at time t, i.e.:

Gaussian case
We first consider the case that the randomness can be modelled as a Gaussian random variable. A similar model is shown in [16].
We choose choose to model the stress as a random noise: where σ 0 and α are constant and ξ is a a Gaussian random variable with zero mean and unit variance. ξ represents the uncertainty related to the value of the stress σ. From Eq (7), the mean value of σ is σ 0 and the variance of σ is equal to α 2 . We assume throughout that the Young's modulus E and the viscosity η are perfectly known. Given that the forcing stress σ for the system is random, the strain is also random. The goal then is to evaluate the derivative of the expected value of the strain with respect to the mean value σ 0 : using the method of MWS. We choose to solve the ODE Eq (5) using an Euler explicit finite difference method with time step δt: Remark. Note that the multiplying term before ξ contains ffiffiffiffi ffi dt p and not δt. This is a 'conforming' discretisation of the stochastic noise term, resulting in a dependence of the variance of the random parameter on the discretisation size. Informally, taking the limit, we can recover the original ODE Eq (5) as: and: Where V½Á is the the variance. Given that E½ ffiffiffiffi ffi dt p x ¼ 0 and E½ð ffiffiffiffi ffi dt p axÞ 2 ¼ a 2 dt, the numerical method in Eq (9) is consistent in the following sense: For this next part, we adopt the same notation as [16]. We denote the strain of the the system at time t and we denote 0 the strain of the the system at time t + δt. Furthermore, we let P () and P( 0 ) be the probability that the strain of the system is and 0 respectively. The propagator W( ! 0 ) must satisfy: Eq (14) means that the probability that the strain of the system is 0 is the sum (integral) of all the probabilities to be at multiplied by the probability that the system passes from state to 0 during δt. Condition Eq (15) comes from the integration of the first condition over 0 .
To derive the analytical expression of the propagator in Eq (18) we start with the fact that x is Gaussian N * (0, δtα 2 ), hence the probability density function is known and must satisfy the following condition: With an integration by substitution from Eq (9), with: we can then show the expression of the propagator, the probability that the system passes from state to 0 during δt is given by [16]: With the propagator in hand we will now see how it is possible to evaluate the Malliavin derivative with the MWS method. To recap, we denote by J() a quantity of interest of our system and we want to compute the derivative of the mean value of this quantity of interest E½J with respect to a parameter " m, in this case σ 0 when σ = σ 0 + αξ. The form of the Malliavin weights q m can be obtained using the following procedure. First, we know that with dP() = P()d, Eq (1) we can write: and by taking the derivative of Eq (19) [11,16,18]: To define a set of rules for updating q m , we differentiate Eq (14) with respect to m: and we obtain the following rule for updating the Malliavin weight: In the example of random stress with σ = σ 0 + ξ, we obtain: For random Young's modulus E = E 0 + ξ we would have the same expression. In the case of random viscosity η = η 0 + ξ we would obtain following the same logic: With the expression for the Malliavin weight Eq (23) in hand we can now implement an algorithm to calculate the derivative. The procedure is very simple; we take Z samples of the evolutions of the stochastic ODE using the explicit Euler scheme whilst simultaneously evolving the Malliavin weight q m . At teach time step Algorithm 1 describes this procedure in more detail.
The deterministic constants are given to be η = 1, E = 1 and we take a time step of δt = 0.01 for the finite difference scheme. We evaluate by the MWS method the derivative of the expected value of with respect to σ 0 for a loading time t 2 [0, T] with T = 30s. In this example the number of realisations is fixed at Z = 20000. We compare the results with the analytical solution which is: We briefly remark that for all numerical results presented in this paper there are two sources of errors committed with respect to the undiscretised problem. The first error is due to the deterministic approximation of the PDE (finite difference or finite element method), and the second due to the stochastic approximation (Monte Carlo estimator). In all cases we drive the error in the deterministic approximation of the PDE far lower than that in the stochastic approximation, such that the error is dominated by the number of realisations Z used in the Monte Carlo estimator.
In Fig 2 we can see that the MWS method gives a good estimation of the Malliavin derivative, particularly in the non-steady state regime t 5s. The relative error and so the global statistical error can become very high when the system reaches a steady state because the value of the sensitivity derivative is constant but the statistical error is compounded after each time step. To address this issue a technique can be employed based on the correlation function [16]. The reader is referred to [16] for further details. We have implemented this correlation correction, which we denote MWS-steady-state, and we can see in Fig 2 that the error is greatly reduced in the steady-state regime. For the numerical examples presented in the following sections, we will consider only the systems undergoing transition or purely steady state systems. Therefore we will not use the MWS-steady-state method again in this paper.

Non-Gaussian case
In this section we explain how to calculate the derivative for non-Gaussian stochastic parameters using the MWS method. The procedure is similar to that shown in the previous part but the rule for updating the Malliavin weight must be modified.
We begin as before by considering uncertainty in the stress σ: where ξ a random variable with probability density function f(x) and c and σ 1 are two constants. We have written Eq (26) in the form of a constant s 0 ¼ s 1 þ cE½z plus a random variable cðz À E½zÞ with zero mean. Then it follows that σ 0 is the mean of the uncertain stress σ. We will use the MWS method to evaluate the derivative of the expected value of the quantity of interest with respect to σ 0 .
.initialisation for z = 0 to Z − 1 do t = 0; .time (t) = 0; .initial condition q σ = 0; .MWS weight for i = 1 to n do Draw realisation of random variable ξ i ; To be able to use the MWS method the probability density function on the parameter must satisfy some regularity conditions, see [11,20,21]. Intuitively, the probability density function must be sufficiently "regular" on R which holds for the Gaussian, log-normal, Beta(α > 1, β > 1) and Gamma(k > 1, θ) distributions. However, a uniform distribution between two values a and b can not be considered "regular" because the probability density function is not differentiable at a and b. Instead, we choose to regularised approximation of a uniform distribution using a Beta(1 + e, 1 + e) random variable with e ( 1. Continuing, we again discretise Eq (5) using an explicit Euler method with time step δt: Alternatively, in the case of uncertainty in the Young's modulus, the discretisation can be written: or, for uncertainty related to the viscosity: The probability density function of the beta distribution, for 0 x 1, and shape parameters α, β > 0, is: The beta function B is a normalisation constant to ensure that the total probability integrates to 1. In general we will evaluate and update the Malliavin weight for the parameter m as: Note that in Eq (31), it is important to check that the condition E½q m ðtÞ ¼ 0 is verified. If E½q m ðtÞ 6 ¼ 0, the updated rule must be corrected. An example of performing this correction is given in the next section entitled extension to second derivative. Finally, we note that for the initial condition we always impose q m (t = 0) = 0.
For the uncertain Young's modulus modelled with a beta distribution, we have: For the uncertain stress with beta distribution, we have: For the uncertain viscosity with beta distribution, we have: These results and further calculations are summarised in Table 1. The Malliavin derivatives of the Kelvin-Voigt model with respect to the mean of the three parameters {σ 0 , η 0 , E 0 } modelled as beta (2,2) distributions are shown in Fig 3. The exact solution is computed semi-analytically using standard integration rules. Good agreement between the MWS and semi-analytical solution is observed for E 0 and σ 0 . For the viscosity η 0 the number of Monte Carlo samples is not sufficient to achieve negligible error, but the overall trend is followed.

Extension to second derivative
The MWS method can also be used to compute the second Malliavin derivative of the expected value of a quantity of interest J. If the quantity of interest does not depend explicitly of the random parameter, the expression given in Eq (3) is valid, but the more general form is the following: In Eq (35), when we want to compute the second derivative the term @J @ " m Â Ã does not vanish because in this case J is the first derivative with respect to " m and therefore depends on the parameter " m in general. By applying Eq (35), we can show for example in the case of uncertain Young's modulus that: with the following updating rule: and: The constant C 2 E and C EE allow to ensure that the expected value of the global Malliavin weight ðq EE ðtÞ þ q E ðtÞ 2 À C EE À C 2 E Þ has an expected value equal to zero. In this specific case we have: The precise specification of the constants depends on the distribution. We compute them analytically or by using standard numerical integration techniques found in e.g. Scipy or Maple. In Fig 4, a comparison between the analytical solution and the MWS method is given for the value of the second derivative depending on time of the expected value of with respect to the Young's modulus. For the sake of example, the problem specification is the same as in previous sections, with the exception that the random variable follows a log-normal(μ, σ) distribution with mean equal to 0.5 and standard deviation equal to 0.25 which corresponds to μ = Calculating the Malliavin derivative of some stochastic mechanics problems −0.804 and σ = 0.473. The analytical solutions for the two constants C EE and C 2 E are in this case: As we can see in Fig 4, the MWS method gives a good approximation for the evaluation of the second derivative with Z = 10 7 realisations.

Extension to random process
In this paper we deal with random noise and in the next section we show numerical results of stochastic mechanics problems where models are defined as PDEs. The probability density function of the random variables used in these examples does not depend on time. Similarly to the Kelvin-Voigt model presented before, we study a time dependent problem in a finite dimensional space by splitting the time interval [0, T] into a finite number of increments. Note that it is also possible to take into account the random noise only at the initial time instead of generating random variables at each time step. It would be possible to extend this work to random process, e.g. by using a Wiener process W(t) which verifies in particular (W(t + δt) − W (t)) * N(0, δt). In this case, even for simple ODEs, it is very difficult to obtain analytical solutions because the probability density function of a random process evolves in time. The Malliavin calculus is very well adapted to address these stochastic problems but requires much more advanced mathematical tools as those presented in this paper. In addition, the Malliavin calculus has the advantage and the specificity that it is possible to directly work in the continuum (infinite dimensional space) to evaluate the sensitivity derivatives. We hope that the first and simple approach restricted to random variables presented in this paper may be of interest to the engineering community and encourage them to investigate the benefits that the Malliavin calculus could provide in the context of stochastic PDEs.

PDE examples
We now turn our attention to models that are defined as PDEs. To solve the deterministic evaluations of the PDEs we use the finite element method. We have chosen to use DOLFIN, part of the FEniCS Project to implement the finite element method solvers [22].

Elastic bar with stochastic Young's modulus
The strong form PDE and boundary conditions of the behaviour of a 1-dimensional elastic bar (see Fig 5) are: We take f = 1, L = 1 and a stochastic Young's modulus: with ξ a random variable with beta(2, 2) distribution. The forward model is described by the following weak residual formulation, find u 2 H 1 D ðO s Þ such that: where the space H 1 D ðO s Þ is the usual Sobolev space of square-integrable functions with squareintegrable weak derivatives on the domain O s ≔ [0, 1] that satisfy the Dirichlet boundary condition u(0) = 0 and H 1 0 ðO s Þ vanish on the whole boundary. We solve the forward model using a piecewise linear finite element method with 1024 cells in the mesh.
The quantity of interest is: The derivative of the expected value of J with respect to the mean value of the Young's modulus " E can be computed analytically in this case: This problem is a stationary (not time-dependent), in contrast to the Kelvin-Voigt model considered previously. However, this stationary problem can be solved using the same techniques. We introduce the concept of pseudo-time, where the system evolves from its initial state at t = 0 to the final solution at time t = T through the a single solution of the PDE Eq (44). Therefore in algorithm 1 we take the pseudo-time step as δt = T and hence n = 1. As before, the Malliavin weight still has initial condition q m (0) = 0.
Finally, the relative error between the MWS method with Z = 5 × 10 5 realisations and the analytical solution is 3.0 × 10 −3 .

Buckling of a hyperelastic beam with stochastic Young's modulus
We study the deformation of a 2D geometrically non-linear hyperelastic beam with stochastic Young's modulus E. We have deliberately designed this problem so that for some values of E the beam undergoes buckling, and for others not.
Consider a hyperelastic body that in its undeformed state occupies the domain O 0 ¼ ½0; L Â ½0; e & R 2 with L = 0.2m and e = 0.01m (see Fig 6), and in its deformed state occupies some (unknown) domain O & R 2 . φ is the map between the material points X in the undeformed domain O 0 and points x in the deformed domain O: The deformation gradient can be written FðXÞ≔ @φ @X . The right Cauchy-Green tensor is then defined as C ≔ F T F.
The Neo-Hookean stored energy density of the body is then: where I 3 ≔ det ðFÞ and and I C ¼ tr ðCÞ. λ and μ are the Lame parameters and can be expressed as a function of the Young's modulus E and Poisson's ratio ν as: We choose to model the Young's modulus as a log-normal random variable with mean value 11 MPa and standard deviation 2 MPa. We take Poisson's ratio as a fixed constant ν = 0.3. Defining the displacement field as u ≔ φ − X and a linear functional f that encodes the external tractions we can characterise the elastic equilibrium displacement field u Ã as the solution to the following minimisation problem: where ½H 1 D ðO 0 Þ 2 is the usual vector-valued Sobolev space of square integrable functions with square integrable derivatives that satisfies the given Dirichlet boundary conditions and dx 0 is a measure on O 0 . We fix the left hand of the beam, u(0, y) = 0 and apply a surface traction in the negative x direction on the right hand of the beam of magnitude f. For one Monte Carlo realisation we solve the non-linear problem using a Newton method from SNES [23] with continuation in the loading parameter f and a third-order backtracking line search. We let the symbolic differentiation capabilities of UFL derive the residual and Jacobian of the forward model for use in the Newton method. We solve the linear systems arising from the Newton iterations using a conjugate gradient method preconditioned using algebraic multigrid (Hypre BoomerAMG [24]) interfaced from PETSc [23].
The quantity of interest is defined as: The Malliavin derivative of E½J with respect to the mean Young's modulus obtained with the MWS method with Z = 3 × 10 3 realisations is: No analytical solution exists for comparison. If we use dolfin-adjoint [25], we can compute the classical derivative of J with respect to the Young's modulus around the mean parameter: In this example the difference between the classical derivative and the Malliavin derivative is quite pronounced. This difference is caused by the presence of an instability (buckling). This instability is not activated when E = E 0 , hence, the classical derivative tells us that J is relatively insensitive to perturbations in the Young's modulus about E 0 . However, the Malliavin derivative tells us that E½J is in fact quite sensitive to changes in the mean of the Young's modulus E 0 . The Malliavin derivative gives us quite a different perspective on the sensitivity of this problem than the classical one.

Incompressible Navier-Stokes equations with stochastic viscosity
We consider the incompressible Navier-Stokes equations on a domain O in R 2 consisting of a pair of momentum and continuity equations: In Eq (54), u refers to the unknown velocity of the fluid, ν is the viscosity of the fluid, p the unknown pressure and f is a given source. The mesh, geometry and boundary conditions for the incompressible Navier-Stokes problem are shown in Fig 7. The viscosity is modelled as a random variable: n ¼ 0:015 þ 0:01ðx À 0:005Þ; ð55Þ with ξ a log-normal distribution with mean equal to 0.5 and standard deviation equal to 0.25. We solve the PDE for a given parameter ν with FEniCS [26] (FE approximation) using Chorin's method with time step δt = 0.01 for t 2 [0, 1], see [27] for more details on the implementation. For one realisation of the viscosity the velocity at time t = 1 s is show in Fig 8. The quantity of interest is the total volume of fluid that exits the right end of the domain: where S p = 0 is the surface with normal vector n on the right side where the pressure is imposed to zero. The derivative of E½J with respect to η 0 obtained with the MWS method for Z = 4 × 10 5 realisations is: @E½J @Z 0 % À 7:9 s=m 2 : ð57Þ No analytical solution exists for comparison. If we use dolfin-adjoint [25], we can compute the derivative of J with respect to the viscosity around the mean parameter:

@J @Z
Z¼Z 0 % À 7:85 s=m 2 : The two sensitivity derivatives are close. In this example, contrary to the hyperelastic example, the Malliavin approach does not give us a particularly different interpretation of the sensitivity.

Conclusion
In this paper we have shown how to calculate the Malliavin derivative using the method of Malliavin Weight Sampling. We have applied the method to some typical mechanics models that can be described by ODEs and PDEs, and solved those models using finite difference and finite element methods. In addition, we have extended the existing practical literature on MWS to non-Gaussian random variables and the calculation of second-order derivatives. We are currently investigating the extension of this work from random parameters to problems with variables modelled as random fields. We are also exploring the use of the Malliavin Calculating the Malliavin derivative of some stochastic mechanics problems derivative in derivative-driven variance reduction methods e.g. [28]. Code showing the calculation of the Malliavin derivative for the examples in this paper is freely available at: https://dx. doi.org/10.6084/m9.figshare.5432722 [29].