A numerical approach to Kolmogorov equation in high dimension based on Gaussian analysis

TL;DR

This paper introduces a novel numerical method for solving high-dimensional Kolmogorov equations by leveraging Gaussian analysis, providing an alternative to Monte Carlo methods with proven convergence.

Contribution

The paper develops a Gaussian-based series expansion approach for high-dimensional Kolmogorov equations, inspired by SPDE structures, with convergence proof and numerical validation.

Findings

Convergent series expansion for solutions in high dimensions.

Numerical tests demonstrate efficiency over traditional methods.

Applicable to infinite-dimensional stochastic systems.

Abstract

For Kolmogorov equations associated to finite dimensional stochastic differential equations (SDEs) in high dimension, a numerical method alternative to Monte Carlo simulations is proposed. The structure of the SDE is inspired by stochastic Partial Differential Equations (SPDE) and thus contains an underlying Gaussian process which is the key of the algorithm. A series development of the solution in terms of iterated integrals of the Gaussian process is given, it is proved to converge - also in the infinite dimensional limit - and it is numerically tested in a number of examples.

Tables

Table 1. Table 1 : Model parameters, 𝐞 𝐞 \mathbf{e} stands for the vector with all components identically 1 1 1 .

Parameter	Value	Description
$d$	$10$	dimension of the problem
$y_0$	$2𝐞$	parameter of the nonlinearity $B$, Polynomial case
$x$	$𝐞$	values where the solution $u(t,x)$ is computed
$\sigma$	$1$	noise
$T$	$1$	final time of computation for $u(t,x_0)$
$H$	$1$	threshold for the initial condition $u_0\left(x\right)$

Equations

$$\{\begin{array}{rl}\mathrm{d}{X}_t& =\left(A{X}_t+B\left(X_t\right)\right)\mathrm{d}t+\sigma \sqrt{Q}\mathrm{d}{W}_t,\\ X_0& =x,\end{array}$$

$$\{\begin{array}{rl}{\partial }_{t}u(t,x)& =\frac{{\sigma }^2}{2}\mathrm{Tr}\left(Q{D}^{2}u(t,x)\right)+⟨Ax+B\left(x\right),Du(t,x)⟩,\\ u(0,x)& =u_0(x),\end{array}$$

$$u(t,x)=\mathbb{E}\left[u_0(X_t^x)\right],$$

$$u(t,x)=\mathbb{E}\left[1_{\left\{\Vert x\Vert >R\right\}}\left(X_t^x\right)\right]=\mathbb{P}\left(\Vert X_t^x\Vert >R\right).$$

$$\{\begin{array}{rl}\mathrm{d}{Z}_t& =A{Z}_t\mathrm{d}t+\sqrt{Q}\mathrm{d}{W}_t,\\ Z_0& =0.\end{array}$$

$$u(t,x)=\sum _{n=0}^{\infty }{v}^n(t,x),$$

$$v^0(t,x)=\mathbb{E}\left[u_0\left(e^{tA}x+\sigma Z_t\right)\right]$$

$$v^n(t,x)=$$

$$\displaystyle\ \mathbb{E}\left[u_{0}\big{(}e^{tA}x+\sigma Z_{t}\big{)}\prod_{i=1}^{n}\left\langle\Xi_{\sigma}(r_{i+1}-r_{i})B\big{(}e^{r_{i}A}x+\sigma Z_{r_{i}}\big{)},Z_{r_{i+1}}-e^{(r_{i+1}-r_{i})A}Z_{r_{i}}\right\rangle\right].$$

$$\Vert v^n(t){\Vert }_{\infty }\le \Vert u_0{\Vert }_{\infty }\Vert B{\Vert }_{\infty }^n{C}_{\delta }^n{t}^{n(1-\delta )}\frac{\Gamma (1-\delta {)}^n}{\Gamma (1+n(1-\delta ))},t>0,$$

$$\partial_{t}u(t,x)=\frac{1}{2}{\rm Tr}\big{(}QD^{2}u(t,x)\big{)}+\big{\langle}Ax+B(x),Du(t,x)\big{\rangle},\quad u(0,\cdot)=u_{0}.$$

$$Q_t={\int }_0^t{e}^{sA}Q{e}^{s{A}^{\ast }}\mathrm{d}s$$

$$\Vert \Lambda (t){\Vert }_{\mathcal{L}(H)}\le C_{\delta }/t^{\delta },t>0.$$

$${\Xi }_{\sigma }(t)=\sigma Q_t^{-1/2}\Lambda (t)=\sigma Q_t^{-1}{e}^{tA};$$

$$Ax=\Delta x,x\in D(A)=H^2(\mathcal{O})\cap H_0^1(\mathcal{O}),$$

$$A{e}_k=-|k{|}^2{e}_k,k\in {\mathbb{N}}^d,$$

$$e_k(\xi )=(2/\pi {)}^{d/2}\sin (k_1{\xi }_1)\cdots \sin (k_d{\xi }_d),\xi \in [0,\pi {]}^d.$$

$$Qx=\sum _{k\in {\mathbb{N}}^d}|k{|}^{-2\alpha }⟨x,e_k⟩e_k,x\in H.$$

$${\rm Tr}(Q_{t})=\sum_{k\in\mathbb{N}^{d}}\frac{1}{2|k|^{2+2\alpha}}\Big{(}1-e^{-2t|k|^{2}}\Big{)}<\infty.$$

$$\Lambda (t)x=\sum _{k\in {\mathbb{N}}^d}\frac{\sqrt{2}|k{|}^{1+\alpha }}{\sqrt{e^{2t|k{|}^2}-1}}⟨x,e_k⟩e_k,x\in H.$$

$$\Vert \Lambda (t){\Vert }_{\mathcal{L}(H)}\le \frac{\sqrt{2C_{\alpha }}}{t^{(1+\alpha )/2}},$$

$$C_{\alpha }=\underset{\theta >0}{\sup }\frac{{\theta }^{1+\alpha }}{e^{2\theta }-1}<+\infty .$$

$$W_t=\sum _{k=1}^{\infty }{W}_t^k{e}_k,t\ge 0,$$

$$\mathrm{d}{Z}_t^x=A{Z}_t^x\mathrm{d}t+\sqrt{Q}\mathrm{d}{W}_t,Z_0^x=x\in H$$

$$Z_t^x=e^{tA}x+W_A(t),t>0,$$

$$W_A(t)={\int }_0^t{e}^{(t-s)A}\sqrt{Q}\mathrm{d}{W}_s.$$

$$\mathbb{E}\big{\langle}h,Q_{t}^{-1/2}W_{A}(t)\big{>}^{2}=|h|_{H}^{2}.$$

$$S_{t}f(x):=\mathbb{E}f(Z_{t}^{x})=\int_{H}f(y)\,N_{e^{tA}x,Q_{t}}({\rm d}y)=\int_{H}f\big{(}e^{tA}x+y\big{)}\,N_{Q_{t}}({\rm d}y).$$

$$\langle h,DS_{t}f(x)\rangle=\mathbb{E}\big{[}f(Z_{t}^{x})\big{\langle}\Lambda(t)h,Q_{t}^{-1/2}\big{(}Z_{t}^{x}-e^{tA}x\big{)}\big{\rangle}\big{]}.$$

$$\Vert D{S}_{t}f{\Vert }_{\infty }\le \Vert f{\Vert }_{\infty }\Vert \Lambda (t){\Vert }_{\mathcal{L}(H)}.$$

$$u(t,x)=(S_{t}u_{0})(x)+\int_{0}^{t}\big{(}S_{t-s}\langle B,Du(s)\rangle\big{)}(x)\,{\rm d}s.$$

Full text

Keywords: Kolmogorov equation, numerical solution, iteration scheme, Gaussian process

1 Introduction

Kolmogorov equations are parabolic equations with a structure directly related to stochastic differential equations (SDEs). The SDEs considered here are in a finite dimensional space but they are inspired by the spatial discretization of stochastic Partial Differential Equations (SPDE). When the noise is additive and the nonlinearity is time-independent, a general form of such SDEs is

[TABLE]

where $x\in\mathbb{R}^{d}$, $\left(W_{t}\right)_{t\geq 0}$ is a Brownian motion in $\mathbb{R}^{d}$ (namely $W_{t}=\left(W_{t}^{1},\ldots,W_{t}^{d}\right)$ where the $W_{t}^{i}$’s are independent real valued Brownian motions), defined on a probability space $\left(\Omega,\mathcal{F},\mathbb{P}\right)$ with a filtration $\left(\mathcal{F}_{t}\right)_{t\geq 0}$, $\sigma$ is a positive real number measuring the strength of the noise, $Q$ is a $d\times d$ positive definite symmetric matrix (the so called covariance matrix of the noise) describing the spatial structure of the noise and $\sqrt{Q}$ is its square root, $A$ is a $d\times d$ matrix and $B:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a function with the degree of regularity specified below. Obviously we could include the scalar $\sigma^{2}$ inside the matrix $Q$ but for certain practical arguments it is useful to distinguish between them. The solution $X_{t}$ is a continuous adapted process in $\mathbb{R}^{d}$. The associated Kolmogorov equation is

[TABLE]

where $u:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, $Du(t,x)$ and $D^{2}u(t,x)$ denote respectively the vector of first partial derivatives and the matrix of second partial derivatives, ${\rm Tr}\left(QD^{2}u(t,x)\right)$ is the trace of the $d\times d$ matrix $QD^{2}u(t,x)$ and $\langle\cdot,\cdot\rangle$ denotes the scalar product in $\mathbb{R}^{d}$. Both for the SDE and the Kolmogorov equation we have used notations which may be adapted to the infinite dimensional case, when $\mathbb{R}^{d}$ is replaced by a Hilbert space (see Section 2 for the general theory); however, the aim of this work is numerical and all objects in the introduction will belong to $\mathbb{R}^{d}$. The link between the Kolmogorov equation and the SDE is

[TABLE]

where $\mathbb{E}$ denotes the mathematical expectation on $\left(\Omega,\mathcal{F},\mathbb{P}\right)$ and $X_{t}^{x}$ is the solution of the SDE above, where the initial condition $x$ is explicitly indicated. Several elements of theory both in finite and infinite dimensions for SDEs and associated Kolmogorov equations can be found in many books, like [6, 8, 9, 16, 17].

Solving the Kolmogorov equation with suitable initial condition $u_{0}$ is a way to compute relevant expected values and probabilities associated to the solution of an SDE. For instance, when $u_{0}(x)=1_{\left\{\left\|x\right\|>R\right\}}$, $u(t,x)$ is the probability that the solution exceeds a threshold $R$:

[TABLE]

The classical method of computing these expected values is the Monte Carlo method (with important variants, see for instance [13, 19]): several realizations of the process $X_{t}^{x}$ are simulated by solving the SDE – typically by Euler method – and then the corresponding values of $u_{0}(X_{t}^{x})$ are averaged. Going beyond this strategy is a fundamental issue, due to its limitations in relevant applications like Geophysics and Climate change projections [15], especially concerning extreme events. The question is whether Kolmogorov equation can be efficiently solved numerically without using the simulation of the SDE. But the problem is that the dimension $d$ is extremely high in these examples and common numerical methods for solution of parabolic equations already require strong computational power when $d=3$, [5, 18]. A grid of $N$ points in $\mathbb{R}$, repeated for all dimensions, give rise to $N^{d}$ grid points, numerically impossible when, for instance, $N=10,\,d=10$ (which still would be an extremely poor approximation). Spectral methods seem to meet the same restrictions: $N^{d}$ is the cardinality of basis elements obtained by tensorization of $N$ basis elements for each space variable.

The problem of dimensionality, the limitations of present methodologies and several motivations are recalled in two recent works [2, 14] which also aim to go beyond Monte Carlo and propose a method based on deep artificial neural networks. We address to these brilliant works for other comments on the problem, see also [7, Introduction]. The approach developed here is however completely different.

Our aim is to take advantage of the probabilistic structure of the problem to devise numerical schemes for the Kolmogorov equation, in particular using Gaussian analysis. We implement a perturbative scheme which links the solution of Kolmogorov equation to a Gaussian process, the solution $Z_{t}$ of the linear stochastic equation

[TABLE]

The idea comes from the theoretical investigations of infinite dimensional Kolmogorov equations associated to SPDEs, see for instance [9, 10]. We modify and adapt that idea giving an explicit formula in terms of a series of Gaussian integrals. We provide here a first glance at the strategy by writing the final formula:

[TABLE]

where

[TABLE]

and for $n\geq 1$

[TABLE]

The matrix $\Xi_{\sigma}(t)$ will be defined in the next sections, see (2.3); it is easily computed by $A$ and $Q$, and it depends on the parameters $t$ and $\sigma$. A theoretical analysis of this series is made, proving the following result.

Theorem 1.1**.**

Assume that $u_{0}$ and $B$ are bounded. Then, under suitable conditions on $A$ and $Q$ (see Hypothesis 2.1 for details), we have the following uniform estimate:

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function, $C_{\delta}>0$ is a constant and $\delta\in(0,1)$ the parameter in (iv) of Hypothesis 2.1.

This theorem sustains the numerical method and stresses the independence on the dimension of certain issues of the method (obviously others, like getting a sample of $Z$, have a cost which increases with $d$). When $\mathbb{R}^{d}$ is replaced by a Hilbert space $H$ (and below we shall formulate the theorem with assumptions in a Hilbert space) it contains also some theoretical novelties with respect to the literature, especially because it provides an explicit formula.

The numerical evaluation of the terms $v^{n}(t,x)$ is made here, in this paper, by Monte Carlo method based on a sample of the process $Z_{t}$ obtained by solving the linear SDE by Euler method. These are the most obvious choices, but other possibilities exist, since $\left(Z_{t}\right)_{t\geq 0}$ is a centered Gaussian process with known covariance function. A main strategy invoked here is to store once for ever a large and accurate sample of $\left(Z_{t}\right)_{t\geq 0}$ (this requires the pair $\left(A,Q\right)$ to be given) and use it later in the formula for different values of the other parameters, $t,x,\sigma,u_{0}$ and even $B$.

This new method is aimed to replace direct Monte Carlo simulations. We should therefore accurately compare them. If the purpose is to make one single computation, classical Monte Carlo wins: the Gaussian method above still requires Monte Carlo simulations of the linear problem, which is less expensive than the nonlinear one but then one has to compute possibly several terms $v^{n}(t,x)$; some experiments clearly show that classical Monte Carlo is less expensive for a comparable degree of precision. The advantage comes when we want to vary parameters, since the Gaussian method for given $\left(A,Q\right)$ allows to store a possibly expensive sample of the process $Z_{t}$ and reuse it for several values of the parameters, just having to compute the averages over the Gaussian sample which give us the terms $v^{n}(t,x)$. On the contrary, classical Monte Carlo method requires to repeat the simulation of the nonlinear problem for each new value of the parameters. By “parameters”, as we have already mentioned above, we mean $t,x,\sigma,u_{0},B$. Let us comment on the interest in changing them.

The interest in changing $t$ is obvious. In certain applications it is necessary to change the initial condition $x$ and compare or collect the results. We have in mind for instance the ensemble methods used in weather prediction where the initial condition is uncertain, a first guess is made on the basis of physical observations, but then the initial condition is perturbed in various directions and the final results averaged by suitable methods. See also [2, 14], where the need to change $\left(t,x\right)$ is stressed.

Changing the strength $\sigma$ of the noise is a very important issue, related also to Large Deviation Theory. We have to advise that the precision of our simulations degenerates as $\sigma\rightarrow 0$, or the number of iterates needed to maintain a reasonable precision blows-up, but at least one can detect some tendency by moving $\sigma$ in a finite range without arriving to too small values.

Concerning the change of function $u_{0}$, unfortunately the main comment is in favor of Monte Carlo: having at disposal a sample of the process $X_{t}^{x}$ immediately gives a way to compute $\mathbb{E}\left[u_{0}(X_{t}^{x})\right]$ for different functions $u_{0}$. Hence the best we can say on this issue is that our formula allows for such computations with a moderate additional effort – but not with an improvement over Monte Carlo.

Finally, changing the nonlinearity $B$ is of theoretical interest for the investigation of the performances of the method, and in applications it may be of interest in those – very common – cases when some parameters of $B$ are not precisely known and different simulations may be useful for comparison or for ensemble averaging methods performed over the range of those parameters.

Let us finally come to a brief description of numerical results. In Section 3, we present some numerical results based on the method proposed here in the finite dimensional settings with $d\geq 10$. The results, even if not fully satisfactory yet, should be compared with the fact that the innovative attempts to solve the Kolmogorov equation in $d>3$ by direct methods, see [7], are often restricted to dimensions smaller than $10$. Large dimension is therefore a very difficult problem that deserves strong effort for improvement, and some of our results – although not in all examples – are quite promising.

As a final comment, let us explicitly mention that the class of Kolmogorov equations studied here is particular, because of the additive and very non-degenerate noise and because we have treated only relatively mild nonlinearities. We have not considered relevant cases from fluid mechanics which have more severe nonlinearities and activation of more scales; after a few initial tests on dyadic models – we point in particular to the recent models on trees which may be very relevant for turbulence theory, see [1, 3, 4] – it was clear that covering these examples with this approach requires further research and improvements. Extension to multiplicative transport noises [11, 12] is another challenging open question.

2 The iteration scheme for Kolmogorov equations on Hilbert spaces

In this section we work in an infinite dimensional separable Hilbert space $H$ and study the iteration scheme for the Kolmogorov equation:

[TABLE]

Here $A:D(A)\subset H\to H$ is an unbounded linear operator, $Q$ is a nonnegative self-adjoint bounded linear operator on $H$, $B:D(B)\subset H\to H$ is a nonlinear measurable mapping and $u_{0}:H\to\mathbb{R}$ is a real valued measurable function. In this section $Q$ plays the role of $\sigma^{2}Q$ to simplify notation. In the following we write $\mathcal{L}(H,H)$ for the Banach space of bounded linear operators on $H$ with the norm $\|\cdot\|_{\mathcal{L}(H)}$.

Throughout this section we assume the following conditions:

Hypothesis 2.1**.**

• (i)

$A:D(A)\subset H\to H$* is the infinitesimal generator of a strongly continuous semigroup $e^{tA}$.*

• (ii)

$Q$* is a nonnegative self-adjoint operator in $\mathcal{L}(H,H)$ satisfying ${\rm Ker}(Q)=\{0\}$, and for any $t>0$ the linear operator*

[TABLE]

is of trace class.

• (iii)

We have $e^{tA}(H)\subset Q_{t}^{1/2}(H)$ for any $t>0$.

• (iv)

Letting $\Lambda(t)=Q_{t}^{-1/2}e^{tA}$, we assume there exist $\delta\in(0,1)$ and $C_{\delta}>0$ such that

[TABLE]

The assumptions (i)–(iii) are quite standard in the literature, see for instance [8, Hypothesis 2.1 and 2.24]. The operator $\Xi_{\sigma}(t)$ appeared in the introduction has the form

[TABLE]

we remark that, in the setting of the introduction, the operator $Q$ in (2.2) should be replaced by $\sigma^{2}Q$ when computing $Q_{t}$. The following example is taken from [8, Example 2.5] which verifies all the assumptions.

Example 2.2**.**

Let $\mathcal{O}=[0,\pi]^{d}$ with $d\in\mathbb{N}$. We choose $H=L^{2}(\mathcal{O})$, and

[TABLE]

where $\Delta$ is the Laplacian operator with Dirichlet boundary condition. $A$ is a self-adjoint negative operator in $H$, and

[TABLE]

where for $k\in\mathbb{N}^{d}$, $|k|^{2}=k_{1}^{2}+\cdots+k_{d}^{2}$ and

[TABLE]

Choose $Q=(-A)^{-\alpha},\,\alpha\in[0,1)$, so that

[TABLE]

For any $t>0$, if $\alpha>d/2-1$, then

[TABLE]

So (ii) is satisfied.

Next, (iii) can be checked by explicit computations. Moreover,

[TABLE]

From this we deduce that

[TABLE]

where

[TABLE]

Thus (iv) holds with $\delta=(1+\alpha)/2\in[1/2,1)$.

We also need the following technical conditions.

Hypothesis 2.3**.**

The initial datum $u_{0}:H\to\mathbb{R}$ and the nonlinear part $B:H\to H$ in (2.1) are bounded and measurable.

This section is organized as follows. In Subsection 2.1, we recall some basic facts in Gaussian analysis on Hilbert space and give the formula for the first term $v^{1}(t,x)$ of the iteration (2.8). We give in Section 2.2 the details for calculating the second term $v^{2}(t,x)$, which will help us to guess and prove the formula for general terms $v^{n}(t,x)$ in Section 2.3. In the last part, we estimate the uniform norm of $v^{n}(t,x)$ and show the convergence of the iteration scheme. The limit is the unique mild solution of (2.1), see Theorem 2.15.

2.1 Some preparations

Let $W$ be a cylindrical Brownian motion on $H$:

[TABLE]

where $\{e_{k}\}_{k\geq 1}$ is a complete orthonormal basis of $H$ and $\{W^{k}\}_{k\geq 1}$ is a family of independent one dimensional standard Brownian motions defined on some probability space $(\Omega,\mathcal{F},\mathbb{P})$. Under the conditions (i) and (ii) in Hypothesis 2.1, the linear SDE

[TABLE]

has a unique solution with the expression

[TABLE]

where $W_{A}(t)$ is the stochastic convolution:

[TABLE]

For any $t>0$, $W_{A}(t)$ is a centered Gaussian variable on $H$ with covariance operator $Q_{t}$. We denote its law by $N_{Q_{t}}({\rm d}y)$. Accordingly, the law of $Z_{t}^{x}$ is denoted as $N_{e^{tA}x,Q_{t}}({\rm d}y)$. Recall that for any $h\in H$, \big{\langle}h,Q_{t}^{-1/2}W_{A}(t)\big{>} is a centered real Gaussian variable with variance

[TABLE]

We shall write $\mathcal{B}(H)$ for the space of bounded measurable functions on $H$ and $C_{b}^{1}(H)$ the space of Fréchet differentiable functions, bounded with bounded derivatives. When $f\in C_{b}^{1}(H)$, its Fréchet derivative will be denoted by $Df$. For any $f\in\mathcal{B}(H)$ and $t\geq 0$, let

[TABLE]

This defines a Markov semigroup on $H$. We have the following important result which implies $S_{t}$ is strong Feller (see [8, Proposition 2.28] for a proof).

Proposition 2.4**.**

Assume the conditions (i)–(iii) in Hypothesis 2.1. Then for all $f\in\mathcal{B}(H)$ and $t>0$, we have $S_{t}f\in C_{b}^{1}(H)$ and for any $h\in H$,

[TABLE]

Moreover,

[TABLE]

Using the semigroup $S_{t}$, the mild formulation of the Kolmogorov equation (2.1) is

[TABLE]

This suggests us to consider the iterative scheme:

[TABLE]

with $u^{0}(t,x)=(S_{t}u_{0})(x)=\mathbb{E}u_{0}(Z_{t}^{x})$. We define $v^{0}(t,x)=u^{0}(t,x)$ and

[TABLE]

then the new functions satisfy the iteration procedure:

[TABLE]

Before concluding this section, we show how to obtain the first term $v^{1}(t,x)$. Since $u_{0}\in\mathcal{B}(H)$, Proposition 2.4 implies $v^{0}(t)\in C_{b}^{1}(H)$ for any $t>0$, and thus $\langle B,Dv^{0}(t)\rangle\in\mathcal{B}(H)$. Denote by $\mathcal{F}_{t}$ the filtration generated by the cylindrical Brownian motion $W_{t}$.

Lemma 2.5**.**

It holds that

[TABLE]

Proof.

Use the property of conditional expectation:

[TABLE]

where the second step follows from the Markov property. Again by the Markov property,

[TABLE]

where the second step is due to (2.5). Substituting this equality into the previous one we obtain the identity. ∎

The above lemma implies

Corollary 2.6**.**

For any $t>0$ and $x\in H$,

[TABLE]

Moreover,

[TABLE]

and

[TABLE]

Proof.

The formula (2.9) follows directly from Lemma 2.5. Next, by the definition (2.8) of the iteration, for any $s>0$ and $y\in H$,

[TABLE]

where the last inequality follows from (2.6). Therefore,

[TABLE]

which yields the estimate on $\|v^{1}(t)\|_{\infty}$. The inequality (2.10) implies that $k^{0}_{s}\in\mathcal{B}(H)$ for all $s>0$, hence by Proposition 2.4, $S_{t-s}k^{0}_{s}\in C_{b}^{1}(H)$ and

[TABLE]

Finally, by (2.6),

[TABLE]

which, together with (2.10), gives us the last estimate. ∎

2.2 The term $v^{2}(t,x)$

In this part, we compute the second term in the iteration to illustrate the ideas. First we prove

Lemma 2.7**.**

One has

[TABLE]

Proof.

By Corollary 2.6, for any $t>0$, $v^{1}(t)\in C_{b}^{1}(H)$ and

[TABLE]

Recall that (2.10) implies $k^{0}_{s}\in\mathcal{B}(H)$, thus by Proposition 2.4,

[TABLE]

According to the proof of Lemma 2.5, we have

[TABLE]

Note that \big{\langle}\Lambda(t-s)B(x),Q_{t-s}^{-1/2}\big{(}Z^{x}_{t-s}-e^{(t-s)A}x\big{)}\big{\rangle} is $\mathcal{F}_{t-s}$-measurable. Substituting this equality into the one above and using the property of conditional expectation, we obtain the desired result. ∎

Now we are ready to present the expression and estimates for the second iteration.

Proposition 2.8**.**

For any $t>0$ and $x\in H$,

[TABLE]

Furthermore,

[TABLE]

and

[TABLE]

Proof.

By Lemma 2.7, for any $s>0$ and $y\in H$,

[TABLE]

We have

[TABLE]

where the second step follows from the Markov property. Therefore,

[TABLE]

Next, by the definition of $k^{1}_{s}$ and the last inequality in Corollary 2.6,

[TABLE]

This immediately implies

[TABLE]

and we obtain the estimate on $\|v^{2}(t)\|_{\infty}$. Moreover, by Proposition 2.4,

[TABLE]

which, combined with (2.11), gives us the second estimate. ∎

2.3 The general terms $v^{n}(t,x)$

In order to do further iteration, we rewrite the formula in Proposition 2.8 as

[TABLE]

Moreover, denoting by $s_{0}=0$, then we have

[TABLE]

From this we can guess the general formulae.

Theorem 2.9**.**

Let $s_{0}=0$. For any $n\geq 1$,

[TABLE]

Moreover,

[TABLE]

and, letting $s_{n+1}=t$,

[TABLE]

Proof.

We proceed by induction. Indeed, in view of the proofs in Section 2.2, we shall also prove inductively the formula

[TABLE]

where $s_{0}=0$ and $s_{n+1}=t$. The discussions in Sections 2.1 and 2.2 show that the assertions on $v$ hold for $n=1,\,2$, and the above formula of $k$ holds with $n=1$. Now we assume the assertions on $v$ (resp. on $k$) hold for $n$ (resp. for $n-1$), and try to prove them in the next iteration.

By the induction hypotheses, we have $v^{n}(s)\in C_{b}^{1}(H)$ for all $s>0$ and thus, by the definition of the iteration (2.8), $k^{n}_{s}\in\mathcal{B}(H)$ with

[TABLE]

where $s_{n+1}=s$. Proposition 2.4 implies $S_{t-s}k^{n}_{s}\in C_{b}^{1}(H)$ for all $s\in(0,t)$, and from the formula

[TABLE]

we deduce readily the estimates on $\|v^{n+1}(t)\|_{\infty}$ and $\|Dv^{n+1}(t)\|_{\infty}$.

Next we prove the formula for $k^{n}_{t}(x)$ (note that the induction hypothesis gives us the expression of $k^{n-1}_{t}(x)$). We have

[TABLE]

where we used Proposition 2.4 in the last step. By the induction hypothesis,

[TABLE]

where $s_{0}=0$ and $s_{n}=s$. Therefore, by the Markov property,

[TABLE]

Inserting this identity into (2.13) and noticing that \big{\langle}\Lambda(t-s)B(x),Q_{t-s}^{-1/2}(Z^{x}_{t-s}-e^{(t-s)A}x)\big{\rangle} is measurable with respect to $\mathcal{F}_{t-s}$, we obtain

[TABLE]

Renaming $s$ as $s_{n}$ gives us the formula of $k^{n}_{t}(x)$ in the new iteration for all $t>0$ and $x\in H$.

Finally we prove the expression for $v^{n+1}(t,x)$. We have

[TABLE]

Using the formula we have just proved for $k^{n}_{s}(y)$ and the Markov property, we can obtain the expression for $v^{n+1}(t,x)$ in a similar way as above. ∎

We give a slightly different formula which is more appropriate for numerical purpose.

Corollary 2.10**.**

For any $n\geq 1$,

[TABLE]

where $r_{n+1}=t$. Accordingly,

[TABLE]

and, setting $r_{0}=0$,

[TABLE]

Proof.

We change variables as follows:

[TABLE]

The domain of integration becomes

[TABLE]

and $s_{i}-s_{i-1}=r_{n+2-i}-r_{n+1-i},\,1\leq i\leq n$. Therefore, by (2.12),

[TABLE]

In the product, letting $j=n+1-i$, we get the desired formula (2.14). The proofs of the two estimates are similar. ∎

Remark 2.11**.**

Due to the convolution structure (2.8), it seems that (2.14) is not suitable for the induction argument in the proof of Theorem 2.9.

2.4 Convergence of the iteration scheme (2.8)

We need the following technical result, where we use the Gamma function $\Gamma(\alpha)$:

[TABLE]

Lemma 2.12**.**

Assume $\delta\in(0,1)$ and $n\geq 1$. Let $r_{0}=0$ and $r_{n+1}=t$. One has

[TABLE]

and

[TABLE]

Proof.

First we prove

[TABLE]

where $B(\alpha,\beta)$ is the Beta function:

[TABLE]

We proceed by induction. For $n=1$, noting that $r_{2}=t$, we change the variable $\theta=r_{1}/t$ and get

[TABLE]

Therefore the equality holds when $n=1$. Now suppose the equality holds for $n-1$, we prove it for $n$. By the induction hypothesis,

[TABLE]

thus, noticing that $r_{n+1}=t$,

[TABLE]

We have, by changing variable $\theta=r_{n}/t$,

[TABLE]

Substituting this result into the previous one gives us the identity (2.15).

Next, it is well known that

[TABLE]

Therefore,

[TABLE]

Combining this with (2.15) we obtain the desired formula.

The proof of the second identity is similar, by first establishing the identity

[TABLE]

We omit the details here. ∎

As a consequence, we have the following estimates.

Corollary 2.13**.**

Under the Hypotheses 2.1 and 2.3, for any $n\geq 0$ and $t>0$,

[TABLE]

and

[TABLE]

Proof.

The case $n=0$ follows directly from (2.6). Combining Lemma 2.12 and Corollary 2.10, we obtain the general cases. ∎

Now we can prove the existence of limit for the iteration scheme (2.8).

Proposition 2.14**.**

Assume the Hypotheses 2.1 and 2.3. For any $T>0$, the series

[TABLE]

converge uniformly on $[0,T]\times H$. Moreover, for any $t_{0}\in(0,T)$, the series

[TABLE]

converge uniformly on $[t_{0},T]\times H$.

Proof.

We only prove the first assertion; the proof of the second one is similar. By Corollary 2.13 and using the ratio test, it is sufficient to show that

[TABLE]

This follows from elementary calculations. Indeed, setting $\alpha=1-\delta$ for simplicity of notation,

[TABLE]

where $(n\alpha)$ is the decimal part of $n\alpha$. Using the simple inequality $\log(1+x)<x$ for all $x\in(-1,0)$, we have

[TABLE]

Hence,

[TABLE]

Note that the first part on the right hand side tends to $-\infty$ as $n\to\infty$, while the last part is uniformly bounded in $n$, thus we conclude the result. ∎

Thanks to Proposition 2.14, we can define the limit

[TABLE]

moreover, for any $t>0$, one has $u(t)\in C_{b}^{1}(H)$ and

[TABLE]

which holds uniformly on $[t_{0},T]\times H$ for any $0<t_{0}<T$. Finally we can prove the main result.

Theorem 2.15**.**

The limit $u(t,x)$ is the unique solution to the Kolmogorov equation (2.1) in the following sense:

• (a)

for any $T>0$, $u(t,x)$ is uniformly bounded for $(t,x)\in[0,T]\times H$, and $u(t)\in C_{b}^{1}(H)$ for any $t>0$;

• (b)

for any $T>0$, one has $\int_{0}^{T}\|Du(t)\|_{\infty}\,{\rm d}t<\infty$;

• (c)

it satisfies the mild formulation (2.7) for any $t>0$ and $x\in H$.

Proof.

Obviously our limit verifies (a). Next,

[TABLE]

Therefore,

[TABLE]

which shows that (b) is also satisfied. Moreover, for any $t>0$ and $x\in H$,

[TABLE]

This implies the integral in the signs of absolute value makes sense.

It remains to check that $u(t,x)$ verify (2.7). By the iteration scheme (2.8), one has, for any $n>1$,

[TABLE]

The left hand side converges uniformly to $u(t,x)$ on $[0,T]\times H$ for any $T>0$. It suffices to show the uniform convergence of the right hand side. We have

[TABLE]

Similarly to the calculations in (2.17), we can show that the right hand side vanishes as $n$ goes to infinity. Therefore we let $n\to\infty$ on both sides of (2.18) and conclude that $u(t,x)$ satisfies (2.7) uniformly in $(t,x)\in[0,T]\times H$.

Finally we prove the uniqueness of solutions. Suppose $u(t,x)$ and $\tilde{u}(t,x)$ are two solutions to (2.1) with the properties (a)–(c). Then, for any $t>0$ and $x\in H$,

[TABLE]

Therefore,

[TABLE]

Moreover, by Proposition 2.4,

[TABLE]

Hence,

[TABLE]

Under Hypothesis 2.1-(iv), there is some $t_{1}>0$ such that $\|B\|_{\infty}\int_{0}^{t_{1}}\|\Lambda(s)\|_{\mathcal{L}(H)}\,{\rm d}s<1$. Then

[TABLE]

Combining this with (2.19) we see that $u(t,x)=\tilde{u}(t,x)$ for any $(t,x)\in[0,t_{1}]\times H$. Next, by the semigroup property, it is easy to show that, for $t\in(0,t_{1}]$,

[TABLE]

Repeating the above procedure we can prove the uniqueness on the interval $[t_{1},2t_{1}]$ and so on. Thus we complete the proof. ∎

3 Numerical Simulations

In this section we propose some experiment of the iteration scheme (2.8) studied in Section 2 in the finite dimensional setting. We have in mind the framework of Example 2.2, i.e. $A=\Delta$. Since we are in the finite dimensional setting this choice corresponds to take $A\in\mathbb{R}^{d}\otimes\mathbb{R}^{d}$ as the diagonal matrix where $A_{k,k}=-k^{2},\,k=1,\ldots,d$. Moreover we consider the matrix $Q=\sigma^{2}I_{d\times d}$ where $I_{d\times d}$ is the identity matrix over $\mathbb{R}^{d}$, and the parameter $\sigma$ will be specified below (see Table 2 for reference parameters).

We will consider two main classes of examples as a benchmark for our approximation scheme. First, we consider the nonlinear vector field

[TABLE]

i.e. we apply the sine function to all the components. This nonlinearity will be the easier one of our examples since it is close to linear, at least for small values of $x$. We will also consider some variation of the previous example, made by

[TABLE]

where $B_{m}\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ is the skew symmetric matrix

[TABLE]

i.e. the Toeplitz matrix with all one above the diagonal and minus one below. This example is more complex than the previous one. It is significant since it deals with skew symmetric matrices, inducing rotations, which are a first simple step in the direction of fluid dynamics. The vector field $B_{m}x$ is also multiplied by the function $\sin(x)$ in order to make example (3.2) nonlinear. Notice that this last example is not covered by our present theory, since it does not satisfy Hypothesis 2.3. However, even if (3.2) is not bounded, it satisfies a linear growth condition. We hope to improve our theory and the generality of the assumptions in such direction in a future research, and limit ourselves to some numerical experiments for the present work.

Second, we consider the following class of polynomial nonlinearities

[TABLE]

where $\overline{y}\in\mathbb{R}^{d}$ is fixed. Note that this example appeals to the one dimensional case

[TABLE]

for which the dynamical system

[TABLE]

has the singleton $\{\overline{y}\}$ as a stable attractor. The reason behind the example (3.3) is the following: it is close to a polynomial nonlinearity, so that it makes a significant test case; at the same time, the normalization by the factor $\left\|\overline{y}\right\|/(\left\|\overline{y}\right\|+\left\|\overline{y}-x\right\|^{p})$ makes it a bounded operator, so that it fulfills Hypothesis 2.3.

In all the examples above we adopt the following choice of initial condition

[TABLE]

where the parameter $H$ is set to $1$ (see Table 2).

3.1 Approximation schemes

Standard Monte Carlo approach.

Since an explicit solution for Equation (1.2) is not available we will always compare to the solution obtained by means of Monte Carlo simulation of the nonlinear process $X^{x}_{t}$:

[TABLE]

where $N_{s}$ is the number of samples considered, and the processes $X^{x,i}_{t},i=1,\dots,N_{s}$ are independent copies of $X^{x}_{t}$. To compute samples of the process $X^{x,i}_{t}$ we use the Euler-Maruyama scheme with a very fine time step in order to get a good approximation to be used as a comparison. The solution computed by (3.4) will always be referred to in what follows as the reference case.

Numerical iteration scheme.

Under our assumption, since $A$ and $Q$ are diagonal, we can rewrite the equations for the processes $Z^{x}_{t}$ and $Z_{t}$ in a simple way: for $k=1,\ldots,d$,

[TABLE]

and

[TABLE]

We remark that, differently from $Z_{t}$, the process $Z^{x}_{t}$ depends also on the parameter $\sigma$, but we do not explicitly write $Z^{x,\sigma}_{t}$ for ease of notation. Note that the process $Z_{t}$ depends only on the operators $A$. This opens the possibility of computing $Z^{x}_{t}$, and hence also $u(t,x)$, for many values of $x$ without repeating the computations for $Z_{t}$. The same reasoning holds for different values of $\sigma$, see Figure 6. Note also that this strategy cannot be applied to the process $X_{t}^{x}$ since in that case the problem is nonlinear.

Once realizations of the process $Z^{x}_{t}$ are computed, we can proceed with the iteration algorithm (2.8). In order to compute numerically the quantity $v^{n}(t,x)$ appearing in Theorem 2.9 one needs to be able to compute first

[TABLE]

Since $A$ and $Q$ are diagonal and explicit (see the beginning of this section), one has

[TABLE]

[TABLE]

and thus,

[TABLE]

Hence, when integrating expression (3.6), by change of variable we have

[TABLE]

Changing variable provides a significant advantage when performing numerical integration. In fact it is more complex to compute $Z_{t-s}^{x}$ than $\Lambda(t-s)$ (resp. $Q_{t-s}^{1/2}$) since $Z^{x}$ is random and hence we would have been obliged to reverse the time for every sample of the process. On the other hand the matrix $\Lambda(t-s)$ (resp. $Q_{t-s}^{1/2}$) is deterministic so that changing time $s\mapsto t-s$ can be done only once.

Moreover, thanks to Corollary 2.10, it is possible to compute $v^{n}(t,x)$ with a single time integration from the previous step. Introduce

[TABLE]

and notice that, due to Equation (2.14), we have

[TABLE]

Since

[TABLE]

once we have computed $I^{n}$, computing $I^{n+1}$ is a matter of a single integration. This is really crucial because, otherwise, by using the direct expression (2.12) in Theorem 2.9, to compute $v^{n}(t,x)$ one should have done an $n$-dimensional numerical integration, independently on the previous iteration.

Stopping conditions.

Since the numerical scheme is iterative and since an exact solution is not available, we adopt a consecutive-iterations stopping condition. At every step we measure the difference between consecutive iterations and stop when this difference is below a certain threshold $tol$. Specifically we adopt two strategies in different situations: when we compute the entire trajectory of $u(t,x)$ for $t\in[0,T]$, we measure

[TABLE]

and stop the iterations if $err(n)<tol$ (see Figures 1 and 2); when we are interested only in $u(T,x)$ for a fixed $T$, then

[TABLE]

and adopt the same stopping rule (Figures 6 and 7).

The entire procedure can be summarized in the following scheme:

3.2 Examples

Here we collect the results obtained, and all the parameters involved in the simulations. Parameters are divided into two categories: those related to the mathematical problem, and those strictly related to the numerical approximations, see Tables 2 and 2. Those are our reference parameters: we will specify each time any modifications.

In all the figures below, when showing the entire trajectory of the solution $u(t,x)$ for $t\in[0,T]$, we also plot the [math]-th order iteration. This corresponds to the solution of the linear case for (1.2), i.e. the Kolmogorov equation with $B\equiv 0$. This will allow us to compare with the linear case, in order to be sure to have introduced a significant nonlinearity into the problem.

Mixed-time-step strategy.

To perform numerical simulation of SDEs and numerical integration we adopt a mixed-time-step strategy. When we compute the reference solution, through the simulation of the process $X^{x}_{t}$, as well as when computing samples of the linear process $Z_{t}$, we adopt a time step $\Delta t$. On the other hand when we perform numerical integration, to compute successive iterations, we adopt a time discretization parameter $dt\gg\Delta t$, see Table 2. This is due to the fact that, in equation (1.1), as well as in (3.5), a coefficient $-k^{2}$ is present in the $k$-th component of the drift of the equation. This coefficient, and hence the Lipschitz constant of the drift, is growing as the square of the dimension $d$ of the problem. This is caused by the intrinsic exponential decay of equation (3.5), which require a high level of precision in computation. Differently, in equation (3.7), part of this exponential decay is absorbed by the convolutional structure of the integration. The limits and what is the proper ratio between $\Delta t$ and $dt$ is a difficult topic. A more precise investigation is needed: for the present paper we only highlight the numerical result obtained, and hope to improve the theoretical counterpart in a future work.

Positive results.

For the simpler test case, the sine case (3.1), see Figure 1, convergence is obtained in five iterations. This is due to the simplicity of the example, as $\sin(x)$ is almost linear near the origin. The situation is different when dealing with some more concrete examples like the polynomial case. In Figure 2, where we use formula (3.3) with $p=2$, we see that the number of iterations to convergence is much bigger (26 in our example). At the same time the difference between the last iteration and the reference case is quite small, comparable to the sine case. However, we notice that the oscillation of the solution computed via our iteration scheme, related to the variance of the estimator, is a bit bigger than that of the reference case. This discrepancy is not completely clear yet, even if we expect it to be due to the low number of samples used to compute averages. In Figure 2 we also add a moving-average smoothing of the solution, to make more perceivable this last intuition.

The same behavior is obtained in the variations of the previous examples. In Figure 3 we see that the same fast convergence as in the sine case, is obtained also in the sine times skew-symmetric case (3.2). The Polynomial cubic case (3.3) with $p=3$ has the same level of complexity as the case with $p=2$, even if it requires a higher number of iterations to obtain convergence, and presents the same type of oscillations.

We also perform the same tests in much higher dimension. In Figure 4 we show the results of the same examples, performed in dimension $d=50$ with $N_{s}=10^{4}$ samples. We see that the number of iterations required to convergence are comparable with result in $d=10$: this confirms the estimate of Corollary 2.13 which is in the infinite dimensional framework and hence is independent of any dimension. The small variations in the number of iterations, as well as the slight increase of the oscillations in the quadratic case, can be explained by the reduction in the number of samples used to compute empirical averages. It is also important to remark that, in the current example, the estimate of Corollary 2.13 is still too rough: by computing the right-hand side of (2.16) one finds that the number $n$ of iterations needed to have $|v^{n}(t,x)|<tol$ is far bigger than what we find in the numerical test (in fact it should be bigger than one hundred).

In Figure 6 we followed a different approach: we fix the test case as the sine times skew-symmetric matrix (3.2), and analyze what is the limit of $u(1,x)$ as $\sigma$ goes to zero. Also in this case the solution computed through the iteration scheme is quite close to the reference case. At the same time, on the right side of Figure 6, we can appreciate the great advantage in time-saving of the iteration scheme. We remark that the plot on the right side is cumulative, meaning that it takes into account the time spent to compute the solution multiple times. In particular, we note that the reference case is a straight line, since the computational time does not depend on the different values of $\sigma$. On the other hand, for the iteration scheme there is a change in the number of iterations for different values of $\sigma$ that justifies the nonlinear shape. Moreover we see that, even including the time of computing samples of the process $Z_{t}$ that can be done only once (since $\sigma$ does not appear in (3.5)), we still have a great advantage in time.

As remarked in the introduction, this kind of advantage is a main feature of the new method proposed here and applies also to the variation of other parameters than $\sigma$. In particular, it applies to the change of initial conditions $x$, one of the most fundamental problems in weather and climate prediction, related to the ensemble forecasting method, see [15, Chapter 6]. Again, Monte Carlo pays linearly with the number of variations of $x$, while our method pays the bulk (i.e. $Z_{t}$ in (3.5)) only once and then (here for the initial conditions) roughly linearly in the number of different $x$’s, but with a linear slope much smaller than the one of Monte Carlo, similarly to the initial slope of Figure 6 right side. We illustrate the interest in varying $x$ by Figure 7 right side, where it is illustrated the relative importance of different variations.

Difficulties with small $\sigma$ and high dimension.

However, not every situation is well behaved as those presented above: in Figure 7 left side, we present the plot for different values of $\sigma$ in the polynomial quadratic case. Here the approximation tends to degenerate for smaller values of $\sigma$ (already around $0.5$). This is due to the higher level of nonlinearity of the polynomial case with respect to (3.2). It is also important to mention that the number of iterations to convergence is really important for what concerns the computational time. In the polynomial quadratic (and also cubic) case, since the number of iterations to convergence is much higher than in the simpler case, the advantage in the computational time is less relevant. Still for what concerns negative results we also show in Figure 5 that, when the dimension grows (left $d=20$, right $d=50$), the sine times skew-symmetric case (3.2) tends to degenerate. Iterations are still converging but the limit is far form the reference solution. This is definitively the most difficult of our examples since it is the only one which mixes strongly all the components and produces a strong energy flux between them. However we also remark that at present time this case is not covered by our theory, but is still relevant since it has the rotational behavior which appeals to fluid dynamics.