Optimal approximation of stochastic integrals in analytic noise model
Andrzej Kałuża Paweł M. Morkisz Paweł Przybyłowicz


Optimal approximation of stochastic integrals in analytic  noise model]
Optimal approximation of stochastic integrals in analytic noise model





A. Kaluza]Andrzej Kałuża
AGH University of Science and Technology,
Faculty of Applied Mathematics,
Al. A. Mickiewicza 30, 30-059 Kraków, Poland
akaluza@agh.edu.pl


P. Morkisz]Paweł M. Morkisz
AGH University of Science and Technology,
Faculty of Applied Mathematics,
Al. A. Mickiewicza 30, 30-059 Kraków, Poland
morkiszp@agh.edu.pl

P. Przybyłowicz]Paweł Przybyłowicz
AGH University of Science and Technology,
Faculty of Applied Mathematics,
Al. A. Mickiewicza 30, 30-059 Kraków, Poland
pprzybyl@agh.edu.pl, corresponding author



We study  approximate stochastic Itô integration of processes belonging to a class of progressively measurable stochastic processes that are Hölder continuous in the rth mean. 

Inspired by increasingly popularity of computations with low precision (used on Graphics Processing Units – GPUs and standard Computer Processing Units – CPU for significant speedup), we introduce a suitable analytic noise model of standard noisy information about X and W. In this model we show that the upper bounds on the  error of the Riemann-Maruyama quadrature are proportional to n^-ϱ+δ_1+δ_2, where n is a number of noisy evaluations of X and W, ϱ∈ (0,1] is a Hölder exponent of X, and δ_1,δ_2≥ 0 are precision parameters for values of X and W, respectively. Moreover, we show that the error of any algorithm based on at most n noisy evaluations of X and W is at least  C(n^-ϱ+δ_1).  Finally, we report numerical experiments performed on both CPU and GPU, that confirm our theoretical findings,  together with some computational performance comparison between those two architectures.
Key words: Wiener process, noisy information,  analytic noise model, optimal approximation, minimal error, GPU
Mathematics Subject Classification: 68Q25,  65C30.


[
    [
    
=====









§ INTRODUCTION

In this paper we investigate the problem of optimal approximation of stochastic integrals of the following form

    ℐ(X,W)=∫_0^T X(t) W(t),

where T>0 and W={W(t)}_t≥ 0 is a one-dimensional Wiener process on some probability space 
(Ω,Σ,ℙ), and we consider integrands X={X(t)}_t∈ [0,T] from  a class of progressively measurable stochastic processes that are Hölder continuous in the rth mean. Such quadrature problems arise, for example, in the context of numerical solutions of stochastic differential equations, see Section 4.4 in <cit.> and <cit.>. Since the exact values of such stochastic integrals are known only in limited cases, an efficient approximation of ℐ(X,W) with the  error as small as possible is of interest. We are aiming at methods that are based only on discrete values of X and W which are, additionally, corrupted with some noise. 

The problem of optimal approximation of stochastic Itô integrals under exact information about X and W has been well studied in the literature, see, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>. Less explored is the approximation of stochastic integrals in the case when  values of X and W  are  corrupted with some noise. The noise may arise from measurement errors, previous computations or simply the floating point numbers representation errors. There is a trend observed in deep learning, where the computations are mostly conducted in lower precision, i.e. using not only single but also half precision. Exemplary performance for NVIDIA V100 graphic card is about 7 TFLOPS for double precision and 14 TFLOPS for single precision. For deep learning purposes, using the nature of the operations and also further lowering the precision to half precision in chosen operations, enabled obtaining up to 112 TFLOPS of performance <cit.>. Hence, it is a huge motivation for analyzing how lowering the precision for other than deep learning applications, influence the computed result.

There are many results on solving problems under noisy information, including the problems such as integrating or approximating  regular functions (see, e.g., <cit.>), L_p approximation of piecewise regular functions (see <cit.>), approximate solving of IVPs (see <cit.>) or PDEs (see <cit.>). 

In this paper we study noisy information for stochastic Itô integration. According to our best knowledge this is the first paper that deals with noisy information for stochastic processes and its application to numerical computation of Itô integrals. In this sense this paper can be seen as the extension of the model presented in <cit.> in the context of SDEs. However, in that paper only drift and diffusion coefficients were corrupted with noise, information about the Wiener process was exact.

We use the Information–Based Complexity framework (see <cit.> and <cit.>). We assume that the algorithms may use only  noisy standard information about the integrand X and the Wiener process W. Namely, let  δ_1, δ_2 ≥ 0 be the precision levels corresponding to the processes X and W, respectively. (The case of δ_1 = δ_2 = 0 corresponds to the exact information.) Available standard information about each coefficient consists of  noisy evaluations of X and W of points t_i, z_i ∈ [0,T]. This means that, for example, for the Wiener process and for a given z_i∈[0,T] an evaluation returns  W̃(z_i) such that | W(z_i) - W̃(z_i)| ≤δ_2(1+|W(z_i)|^s) for some s≥ 0. In the context of computations performed on GPUs, this can be interpreted as the standard relative error. (Detailed description of the noisy  information is given in Section 2.) From the reasons that become clear in Section 2 we refer to the model presented as to analytical noise model.

The error of an algorithm is measured in the r-th mean  maximized over the class of input data X and over all permissible information about (X,W) with the given precisions δ_1, δ_2≥ 0. In the model we show that the upper bounds on the  error of the Riemann-Maruyama quadrature are proportional to n^-ϱ+δ_1+δ_2, where n is a number of noisy evaluations of X and W, ϱ∈ (0,1] is a Hölder exponent of X, and δ_1,δ_2≥ 0 are precision parameters for values of X and W, respectively. Moreover, we show that the error of any algorithm based on at most n noisy evaluations of X and W is at least  C(n^-ϱ+δ_1). We also present numerical experiments that confirm our theoretical findings. As we perform the similar operations for multiple trajectories of considered processes, the proposed algorithm is highly parallel. Hence, the usage of GPU for computation acceleration is very natural. There is a vast list of problems where employing GPUs gives significant speedups <cit.>, including matrix operations <cit.>, bioinformatics <cit.> or solving ordinary or random differential equations <cit.>.


The paper is organized as follows. Section 2 contains basic notion and definitions. The Riemann-Maruyama quadrature rule 𝒜_n^RM for the perturbed information and upper estimate on its error (Theorem 1) are presented in Section 3. Section 4 consists of some lower bounds (Lemma <ref>, Proposition <ref>). This leads to the conclusion that the algorithm 𝒜_n^RM is optimal (Theorem <ref>). Section 5 reports
numerical experiments performed for the GPU implementation of the algorithm.








§ PRELIMINARIES

Denote ℕ={1,2,…}. Let  W = {W(t)}_t≥ 0 be the standard, one-dimensional Wiener process defined on a complete probability space  (Ω,Σ, ℙ). By {Σ_t}_t≥ 0 we denote a filtration, satisfying the usual conditions, such that W is a Wiener process with respect to {Σ_t}_t≥ 0.

For a random variable Y:Ω→ℝ we write Y_q=(𝔼|Y|^q)^1/q, q≥ 2. Moreover, by ℒ we mean the following differential operator

    ℒ=∂/∂ t+1/2∂^2/∂ y^2.

For r,q∈ [2,+∞), q≥ r, L≥ 0, ϱ∈ (0,1] we consider the following class F^ϱ,r,q_L of stochastic processes X={X(t)}_t∈ [0,T] 


    F^ϱ,r,q_L={ X:[0,T]×Ω→ℝ |     X  is {Σ_t}_t≥ 0-progressively measurable,
       sup_t∈ [0,T] |X(t)|_q≤ L,
       X(t)-X(z)_r≤ L|t-z|^ϱ,  t,z∈[0,T]}.


Let us recall that by Theorem 33 from Chapter IV in <cit.> for X being {Σ_t}_t≥ 0-progressively measurable, the process {sup_0≤ z≤ t|X(z)|}_t∈ [0,T] is also {Σ_t}_t≥ 0-progressively measurable. Hence, sup_t∈ [0,T] |X(t)|
is (Σ_T-measurable) random variable. Moreover, the processes from the class F^ϱ,r,q_L are Itô integrable, see, for example, <cit.> and <cit.>. 

The numbers r,q,L,ϱ,T   will be called parameters of the class F^ϱ,r,q_L. Except for T the parameters are not known
and the algorithm presented later on will not use them as input parameters.

In order to define suitable model of computation under inexact information about X, W we need to introduce the following auxiliary classes.
Let

    𝒦^1=     {p:[0,T]×ℝ→ℝ |  p(·,·)-Borel measurable,
         |p(t,y)|≤ 1+|y|  for all t∈ [0,T], y∈ℝ}.


For s∈ [0,+∞)  we define

    𝒦^2_s=    {p:[0,T]×ℝ→ℝ |  p∈ C^1,2([0,T]×ℝ),
        max{|∂ p/∂ t(t,y)|,|∂ p/∂ y(t,y)|,|∂^2 p/∂ y^2(t,y)|}≤ 1+|y|^s  for all t∈ [0,T], y∈ℝ},



    𝒦̅^2_s=     {p:[0,T]×ℝ→ℝ |  p∈ C^1,1([0,T]×ℝ),
        max{|∂ p/∂ t(t,y)|,|∂ p/∂ y(t,y)|}≤ 1+|y|^s  for all t∈ [0,T], y∈ℝ},

(if s=0 then we set |y|^s:=0 for all y∈ℝ), 

and for α,β∈ (0,1] 

    𝒦^3_α,β=    {p:[0,T]×ℝ→ℝ |  |p(t,x)-p(z,y)|≤ |t-z|^α+|x-y|^β
         for all t,z∈ [0,T], x,y∈ℝ}

We have that 𝒦^2_0⊂𝒦^3_1,1. In the sequel, the classes above will allow us to to model, at least in some sense, the influence of the regularity of the noise on the error bound.

For δ_1,δ_2≥ 0 we define

    V_X(δ_1)={X̃ |  ∃_p_X∈𝒦^1:∀_(t,ω)∈ [0,T]×Ω X̃(t,ω) = X(t,ω)+δ_1· p_X(t,X(t,ω))},

and the classes of disturbed Wiener process

    𝒲_s(δ_2)={W̃ |  ∃_p_W∈𝒦^2_s:∀_(t,ω)∈ [0,T]×Ω W̃(t,ω) = W(t,ω)+δ_2· p_W(t,W(t,ω))},


    𝒲̅_s(δ_2)={W̃ |  ∃_p_W∈𝒦̅^2_s:∀_(t,ω)∈ [0,T]×Ω W̃(t,ω) = W(t,ω)+δ_2· p_W(t,W(t,ω))},


    𝒲_α,β(δ_2)={W̃ |  ∃_p_W∈𝒦^3_α,β:∀_(t,ω)∈ [0,T]×Ω W̃(t,ω) = W(t,ω)+δ_2· p_W(t,W(t,ω))}.

Since we impose some functional structure on corrupting functions, we refer to this model as to analytic noise model. (See also Remark <ref> for possible alternative approach.)


For X∈ F^ϱ,r,q_L let X̃∈ V_X(δ_1) and W̃∈𝒲(δ_2) where  𝒲∈{𝒲_s,𝒲_α,β}. We assume that the algorithm is based on discrete noisy information about X and W. Hence, a vector of noisy information has the following form

    𝒩(X̃,W̃)=[X̃(t_0),X̃(t_1),…,X̃(t_i_1-1),W̃(z_0),W̃(z_1),…,W̃(z_i_2-1)],

where i_1,i_2∈ℕ. Moreover, t_0,t_1,…,t_i_1-1∈ [0,T] and z_0, z_1,…, z_i_2-1∈ [0, T ] are
given time points. Hence, the information is nonadaptive (see <cit.> and <cit.> for more discussion on adaptive and nonadaptive information). We assume that t_i ≠ t_j, z_i ≠ z_j for all i ≠ j. The total number of (noisy) evaluations of X and W is l=i_1+i_2.

An algorithm 𝒜 using 𝒩(X̃,W̃), that approximates ℐ(X,W), is of the form

    𝒜(X̃,W̃,δ_1,δ_2)=φ(𝒩(X̃,W̃)),

where

    φ:ℝ^i_1+i_2→ℝ,

is a Borel measurable mapping.

For a given n∈ℕ we denote by Φ_n a class of all algorithms of the form (<ref>) for which the total number of evaluations l
is at most n.

For a fixed X∈ F^ϱ,r,q_L the error of 𝒜∈Φ_n is defined as

    e^(r)(𝒜,X,𝒲,δ_1,δ_2)=sup_(X̃,W̃)∈ V_X(δ_1)×𝒲(δ_2)ℐ(X,W)-𝒜(X̃,W̃,δ_1,δ_2)_r,
 
where 𝒲∈{𝒲_s,𝒲̅_s,𝒲_α,β}.





The worst case error of 𝒜 in 𝒢 is given by

    e^(r)(𝒜,𝒢,𝒲,δ_1,δ_2)=sup_X∈𝒢e^(r)(𝒜,X,𝒲,δ_1,δ_2)

where 𝒢 is a subclass of F^ϱ,r,q_L. Finally, the nth minimal error is defined as

    e_n^(r)(𝒢,𝒲,δ_1,δ_2)=inf_𝒜∈Φ_ne^(r)(𝒜,𝒢,𝒲,δ_1,δ_2).

The aim is to develop an optimal algorithm and its efficient implementation by using GPUs.

Unless otherwise stated, all constants appearing in this paper (including those in the '𝒪', 'Ω', and 'Θ' notation) will only
depend on the parameters of the respective classes. Furthermore, the same symbol may be used for different constants.








§ THE RIEMANN-MARUYAMA QUADRATURE FOR NOISY INFORMATION

Let n∈ℕ and

    0=t_0<t_1<…<t_n=T,

be an arbitrary discretization on [0,T]. We denote by Δ t_i=t_i+1-t_i, i=0,1,…,n-1. We define the Riemann-Maruyama quadrature that use noisy evaluations of X and W by

    𝒜^RM_n(X̃,W̃)=∑_i=0^n-1X̃(t_i)· (W̃(t_i+1)-W̃(t_i)),

where (X̃,W̃)∈ V_X(δ_1)×𝒲(δ_2) for 𝒲∈{𝒲_s,𝒲̅_s,𝒲_α,β}.
It is easy to see that the information cost of computing 𝒜^RM_n(X̃,W̃) is 2n noisy
evaluations of X and W. The combinatory cost consists of 𝒪(n) arithmetic operations.

The aim of this section it to prove the following result.
 
 Let us assume that ϱ∈ (0,1] and r≥ 2.


  (i) Let s≥ 0 and q∈ (r,+∞). There exists a positive constant C, depending only on the parameters of the class F^ϱ,r,q_L and s, such that for all n∈ℕ, δ_1,δ_2≥ 0, X∈ F^ϱ,r,q_L, (X̃, W̃)∈ V_X(δ_1)×𝒲_s(δ_2) it holds

    ℐ(X,W)-𝒜^RM_n(X̃,W̃)_r≤ C(max_0≤ i≤ n-1(Δ t_i)^ϱ+δ_1+δ_2+δ_1·δ_2).


  (ii)Let s≥ 0 and q∈ (r,+∞). There exists a positive constant C, depending only on the parameters of the class F^ϱ,r,q_L and s, such that for all n∈ℕ, δ_1,δ_2≥ 0, X∈ F^ϱ,r,q_L, (X̃, W̃)∈ V_X(δ_1)×𝒲̅_s(δ_2) it holds

    ℐ(X,W)-𝒜^RM_n(X̃,W̃)_r   ≤    C(max_0≤ i≤ n-1(Δ t_i)^ϱ+δ_1+
          δ_2(1+δ_1)·(1+∑_i=0^n-1(Δ t_i)^1/2))


  (iii) Let α,β∈ (0,1] and q∈ [r,+∞). There exists a positive constant C, depending only on the parameters of the class F^ϱ,r,q_L and α, β, such that for all n∈ℕ, δ_1,δ_2≥ 0, X∈ F^ϱ,r,q_L, (X̃, W̃)∈ V_X(δ_1)×𝒲_α,β(δ_2) it holds

    ℐ(X,W)-𝒜^RM_n(X̃,W̃)_r    ≤    C(max_0≤ i≤ n-1(Δ t_i)^ϱ+δ_1+
          δ_2(1+δ_1)·∑_i=0^n-1((Δ t_i)^α+(Δ t_i)^β/2)).
  
  

Proof. Let X̃∈ V_X(δ_1). We first show (<ref>), where W̃∈𝒲_s(δ_2). Let the process Z={Z(t)}_t∈ [0,T] be defined as Z(t)=p_W(t,W(t)). Then, by the Itô formula we get that

    Z(t)=M(t)+V(t),    t∈ [0,T],

where


    V(t)=∫_0^tℒp_W(z,W(z)) z,


    M(t)=∫_0^t∂ p_W/∂ y(z,W(z)) W(z).

We stress that {V(t)}_t∈ [0,T] is continuous process with bounded variation, while {M(t)}_t∈ [0,T] is continuous martingale with respect to the filtration {Σ_t}_t≥ 0. Hence, the process Z is continuous semimartingale. 

We denote by 

    Δ Y_i = Y(t_i+1)-Y(t_i)       i=0,1…,n-1,

for Y∈{W,Z} and for all t∈ [0,T]

    X̂_n(t)=∑_i=0^n-1X(t_i)·1_(t_i,t_i+1](t),


    p̂_X,n(t)=∑_i=0^n-1p_X(t_i,X(t_i))·1_(t_i,t_i+1](t).

Note that {X_n(t)}_t∈ [0,T] and {p_X,n(t)}_t∈ [0,T] are {Σ_t}_t≥ 0-progressively measurable simple processes.
Since Z and W are continuous semimartingales, by Property (v) at page 110 in <cit.> we can write the algorithm 𝒜^RM_n as follows

    𝒜^RM_n(X̃,W̃)   =   ∑_i=0^n-1(X(t_i)+δ_1· p_X(t_i,X(t_i)))· (Δ W_i+δ_2·Δ Z_i)
       =   ∫_0^T X̂_n(t) W(t)+δ_1∫_0^T p̂_X,n(t) W(t)
          +δ_2∫_0^TX̂_n(t) Z(t)+δ_1·δ_2∫_0^T p̂_X,n(t) Z(t).

We thus obtain

    ℐ(X,W)-𝒜^RM_n(X̃,W̃)_r≤∑_i=1^4 A_i,n,

where 

    A_1,n=∫_0^T (X̂_n(t)-X(t)) W(t)_r,


    A_2,n=δ_1·∫_0^T p̂_X,n(t) W(t)_r,


    A_3,n=δ_2·∫_0^T X̂_n(t) Z(t)_r,


    A_4,n=δ_1·δ_2·∫_0^T p̂_X,n(t) Z(t)_r.

By the Burkholder inequality and the Hölder continuity of X in rth mean  we get

    A_1,n≤ C_1·(∑_i=0^n-1∫_t_i^t_i+1𝔼|X(t)-X(t_i)|^r  t)^1/r≤ C_2max_0≤ i≤ n-1(Δ t_i)^ϱ,

and

    A_2,n   ≤   δ_1· C_3·(𝔼∫_0^T|p̂_X,n(t)|^r t)^1/r=δ_1· C_3·(∑_i=0^n-1𝔼|p_X(t_i,X(t_i))|^r·Δ t_i)^1/r
       ≤   δ_1· C_4(1+sup_t∈ [0,T]|X(t)|_r)≤δ_1· C_4(1+sup_t∈ [0,T]|X(t)|_q)
       ≤   δ_1· C_5(1+L),

since p_X is of at most linear growth. 
Since Z=M+V, from Definition 5.7 at page 109 in <cit.> we obtain

    A_3,n≤δ_2· (B_1,n+B_2,n),

where

    B_1,n=∫_0^TX̂_n(t)  M(t)_r,


    B_2,n=∫_0^TX̂_n(t)  V(t)_r.

Note that 

    sup_t∈ [0,T]|X̂_n(t)|≤sup_t∈ [0,T]|X(t)|,
 

    sup_t∈ [0,T]|p̂_X,n(t)|≤ 1+sup_t∈ [0,T]|X(t)|,

and

    𝔼(∫_0^T|∂ p_W/∂ y(t,W(t))|^r t)^q/q-r≤ C_6,


    𝔼(∫_0^T|ℒp_W(t,W(t))|^r t)^q/q-r≤ C_7,

since ∂ p_W/∂ y and ℒp_W(t,y) are of at most linear growth. (The constants C_6,C_7  depend only on T, s, r, and q.) Hence, by the associativity property (see, for example, Property (ii) at page 109 in <cit.>), Burkholder and Hölder inequalities we get

    B_1,n   =   ∫_0^TX̂_n(t)·∂ p_W/∂ y(t,W(t))  W(t)_r≤ C_8(𝔼∫_0^T|X̂_n(t)|^r·|∂ p_W/∂ y(t,W(t))|^r t)^1/r
       ≤    C_8sup_t∈ [0,T]|X̂_n(t)|_q·(𝔼(∫_0^T|∂ p_W/∂ y(t,W(t))|^r t)^q/q-r)^1/r-1/q
    
    
       ≤    C_9sup_t∈ [0,T]|X(t)|_q≤ C_9L,

and

    B_2,n   =   ∫_0^TX̂_n(t)·ℒp_W(t,W(t))  t_r≤ C_10(𝔼∫_0^T|X̂_n(t)|^r·|ℒp_W(t,W(t))|^r t)^1/r
    
       ≤    C_11sup_t∈ [0,T]|X̂_n(t)|_q·(𝔼(∫_0^T|ℒp_W(t,W(t))|^r t)^q/q-r)^1/r-1/q
       ≤    C_12sup_t∈ [0,T]|X(t)|_q≤ C_12L,

where C_9,C_12 depend only on T, s, q, r and L. Therefore, by (<ref>), (<ref>), and (<ref>) we arrive at

    A_3,n≤ C_13·δ_2.

By proceeding analogously as for A_3,n we obtain 

    A_4,n≤ C_14·δ_1·δ_2.
 
Combining (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) we get (<ref>), which ends the proof of (<ref>).


We now justify (<ref>) and (<ref>). In this cases the process Z is not necessarily a semimartingale. Hence, we use the following decomposition of 𝒜_n^RM, that follows directly from (<ref>),

    𝒜^RM_n(X̃,W̃)=∫_0^T X̂_n(t) W(t)+δ_1∫_0^T p̂_X,n(t) W(t)
          +δ_2∑_i=0^n-1X(t_i)·Δ Z_i+δ_1·δ_2∑_i=0^n-1p_X(t_i,X(t_i))·Δ Z_i.
 
We have that

    ℐ(X,W)-𝒜^RM_n(X̃,W̃)_2≤∑_i=1^4 D_i,n,

where 

    D_1,n=A_1,n,


    D_2,n=A_2,n,


    D_3,n=δ_2·∑_i=0^n-1X(t_i)·Δ Z_i_r,


    D_4,n=δ_1·δ_2·∑_i=0^n-1p_X(t_i,X(t_i))·Δ Z_i_r.

For D_1,n and D_2,n we use the bounds (<ref>), (<ref>) obtained for A_1,n and A_2,n, respectively. However, for D_3,n and D_4,n we have to differ between the case when W̃∈𝒲̅_s(δ_2) and W̃∈𝒲_α,β(δ_2).

Let W̃∈𝒲̅_s(δ_2). Since in this case p_W∈𝒦̅^2_s we get, by the mean value theorem, that 


    |p_W(t,x)-p_W(z,y)|≤ C_1(|t-z|+(1+|x|^s+|y|^s)· |x-y|),

for all t,z∈ [0,T] and x,y∈ℝ, where C̅>0 depends only on T and s. This implies

    |Δ Z_i|^rγ≤ C_2(Δ t_i)^rγ+C_3(1+|W(t_i)|^rsγ+|W(t_i+1)|^rsγ)·|Δ W_i|^rγ,

for i=0,1,…,n-1, where γ=q/(q-r). Hence, by the Hölder inequality

    X(t_i)·Δ Z_i_r≤sup_t∈ [0,T]|X(t)| _q·(𝔼|Δ Z_i|^rγ)^1/(rγ),

and

    𝔼|Δ Z_i|^rγ   ≤    C_3(𝔼(1+|W(t_i)|^rsγ+|W(t_i+1)|^rsγ)^2)^1/2·(𝔼|Δ W_i|^2rγ)^1/2
          +C_2(Δ t_i)^rγ≤ C_2(Δ t_i)^rγ+C_4(Δ t_i)^rγ/2,

since W has all absolute moments bounded. Therefore,

    X(t_i)·Δ Z_i_r≤ C_5 sup_t∈ [0,T]|X(t)| _q·(Δ t_i+(Δ t_i)^1/2),

which implies that

    D_3,n=δ_2·∑_i=0^n-1X(t_i)·Δ Z_i_r≤δ_2·∑_i=0^n-1X(t_i)·Δ Z_i_r≤ C_6·δ_2·(T+∑_i=0^n-1(Δ t_i)^1/2).

For D_4,n we proceed analogously as for D_3,n.  This gives (<ref>).

Finally, let W̃∈𝒲_α,β(δ_2). Since X(t_i) and Δ W_i are independent, we have in this case that

    D_3,n=δ_2·∑_i=0^n-1X(t_i)·Δ Z_i_r
    
    	   ≤   δ_2·∑_i=0^n-1 |X(t_i)| _r·(Δ t_i)^α+|Δ W_i|^β_r
       ≤   δ_2·sup_t∈ [0,T]|X(t)|_r·∑_i=0^n-1((Δ t_i)^α+|Δ W_i|^β_r)
       ≤   δ_2· L·∑_i=0^n-1((Δ t_i)^α+m^β_rβ(Δ t_i)^β/2),

where m_rβ=Z_rβ and Z is a standard normal random variable with mean zero and variance equal to 1. 




For D_4,n we proceed analogously as for D_3,n. This ends the proof.     ▪

Directly from Theorem <ref> we have the following corollary that states the worst-case error of the algorithm 𝒜^RM_n  in the class F^ϱ,r,q_L.

 Let ϱ∈ (0,1], r≥ 2, and let us consider the Riemann-Maruyama quadrature 𝒜^RM_n based on the equidistant mesh t_i=iT/n, i=0,1,…,n. 

	
  (i) Let s≥ 0 and q∈ (r,+∞). Then
				
    e^(r)(𝒜^RM_n,F^ϱ,r,q_L,𝒲_s,δ_1,δ_2)=𝒪(n^-ϱ+δ_1+δ_2),

				as n→ +∞ and max{δ_1,δ_2}→ 0+.
	
  (ii) Let s≥ 0 and q∈ (r,+∞). Then
				
    e^(r)(𝒜^RM_n,F^ϱ,r,q_L,𝒲̅_s,δ_1,δ_2)=𝒪(n^-ϱ+δ_1+δ_2· (1+n^1/2)),

				as n→ +∞ and max{δ_1,δ_2}→ 0+.			
	
  (iii) Let α,β∈ (0,1] and q∈ [r,+∞). Then
	
    e^(r)(𝒜^RM_n,F^ϱ,r,q_L,𝒲_α,β,δ_1,δ_2)=𝒪(n^-ϱ+δ_1+δ_2· n^1-min{α,β/2}),

	as n→ +∞ and max{δ_1,δ_2}→ 0+.
	


Let us comment on the result obtained so far.
 As we can see from Theorem <ref> and Corollary <ref> domination of the noise term become more on more visible as the regularity of  disturbing functions p_W is decreasing. 



	
	We considered the setting that we called analytic noise model, since we assumed certain form  of the noise via disturbance function p. Of course another approach is possible. Namely, one can assume that the exact values of X are corrupted by noise in the following way
	
    X̃(t_i)=X(t_i)+ε_i,

	where  (ε_i)_i=0,1,…,n are σ(⋃_t≥ 0Σ_t)-measurable random variables. Preliminary estimates indicate that it is possible to achieve upper bounds like in Theorem <ref>, under certain assumptions on the discrete-time process (ε_i)_i=0,1,…,n. We postpone this problem to our future work.








§ LOWER BOUNDS AND OPTIMALITY OF THE RIEMANN-MARUYAMA QUADRATURE
 
In this section we investigate lower bounds on the worst-case error of an arbitrary
algorithm from Φ_n  and, in particular cases, we establish optimality of the Riemann-Maruyama algorithm 𝒜^RM_n. We concentrate on the class 𝒲_s of noisy evaluations of W. Essentially sharp lower bounds in the classes 𝒲̅_s and
 𝒲_α,β are left as an open problem.
 
The following lemma follows directly from (91) in <cit.>, where the lower bound on the error for approximating Itô integrals of deterministic functions from the Hölder class has been established.
 

Let ϱ∈ (0,1], r≥ 2, q∈ (r,+∞), and s≥ 0. Then 

    e^(r)(F^ϱ,r,q_L,𝒲_s,δ_1,δ_2) = Ω( max{n^-ϱ, δ_1}).











From Corollary <ref> and Lemma <ref> we get the  main result of this paper.
 

Let ϱ∈ (0,1], r≥ 2, q∈ (r,+∞), and s≥ 0. Then  the nth minimal error satisfies

    e^(r)_n(F^ϱ,r,q_L,𝒲_s,δ_1,0)=Θ(n^-ϱ+δ_1),

and
	
    e^(r)_n(F^ϱ,r,q_L,𝒲_s,δ_1,δ_1)=Θ(n^-ϱ+δ_1),

	as n→ +∞ and δ_1→ 0+.
An optimal algorithm is the Riemann-Maruyama quadrature 𝒜^RM_n based on the equidistant discretization t_i=iT/n, i=0,1,…,n. 


The results above hold for particular values of the precision parameter δ_2, namely, for δ_2=0 and δ_2=δ_1. In general case preliminary estimates suggest that in order to establish dependence of the lower bounds also on δ_2 completely new technique is required. (We stress that the results from <cit.> are not applicable here, since we consider a different model of noise.) Nevertheless, for the algorithm 𝒜^RM_n we have the following sharp (worst-case) error bounds in the case of arbitrary δ_2.

  
Let ϱ∈ (0,1], r≥ 2, q∈ (r,+∞), s≥ 0, and let us consider the Riemann-Maruyama quadrature 𝒜^RM_n based on the equidistant mesh t_i=iT/n, i=0,1,…,n.
	 Then
				
    e^(r)(𝒜^RM_n,F^ϱ,r,q_L,𝒲_s,δ_1,δ_2)=Θ(n^-ϱ+δ_1+δ_2),

	as n→ +∞ and max{δ_1,δ_2}→ 0+.			

Proof. Upper bounds in (<ref>) follows directly from Corollary <ref>. In the case when δ_1≥ 0=δ_2, the lower bound Ω(max{n^-ϱ,δ_1}) again follows from (91) in <cit.>.

Now we consider the case δ_2≥ 0=δ_1. Let us take

    X_0(t,ω)=L,    (t,ω)∈ [0,T]×Ω,

and take W̃∈𝒲_s(δ_2) of the following form

    W̃(t)=W(t)+δ_2· t.

We get that

    ℐ(X_0,W)=L· W(T),

and

    𝒜^RM_n(X_0,W̃)=L∑_i=0^n-1Δ W_i+Lδ_2∑_i=0^n-1Δ t_i=ℐ(X_0,W)+Lδ_2 T,

which gives 

    ℐ(X_0,W)-𝒜^RM_n(X_0,W̃)_r=LTδ_2.

This implies the thesis.    ▪

In the case of exact information (i.e., δ_1=δ_2=0)  we know, by the results of <cit.>, that even randomized adaptive information does not help, and the rate n^-ϱ is optimal.







§ NUMERICAL RESULTS
 
We present results for the  Riemann-Maruyama quadrature 𝒜_n^RM. There will be four exemplary problems presented, where for the first one we know the exact solution and for the others  we need to assume some convergence of the analyzed method to estimate the obtained error. In the end of this section some practical guides on how to implement the algorithm efficiently using GPUs will be presented together with the discussion about the obtainable speedup of using such architecture.




















































 §.§ Problems

For the test purposes we analyze following integration problem

    ℐ(X,W)=∫_0^T X(t) W(t)

and we consider the following four examples

    X_1(t)    = W(t),    ℐ(X_1,W) = 1/2W^2(T) - 1/2T, 
    
    	X_2(t)     = W_2(t),
    
    	X_3(t)    = f(S(t))S(t),      f(x) = max{0, K-x},      S(t)=S(0)e^-1/2σ^2t+σ W(t)
        K = 9, σ = 1  , S(0) = 1 ,
    
    	X_4(t)     = N(t) e^W(t),

where N={N(t)}_t∈ [0,T] is a Poisson process with insensitivity λ= 5 and W_2={W_2(t)}_t∈ [0,T] is a standard one-dimensional Wiener process, both independent from W. We also apply 𝒜_n^RM to the weak approximation of the following scalar SDE

    {[ dY(t)=μ Y(t) t+W_2(t) W(t),                   t∈ [0,T],;                     Y(0)=0, ].

where μ=3.  The exact solution of (<ref>) leads to the quadrature problem, since

    Y(T)=ℐ(X,W)=∫_0^T e^μ (T-t)W_2(t) W(t),

where X(t)=e^μ (T-t)W_2(t), see <cit.>, <cit.>. We use GPU implementation of the Riemann-Maruyama quadrature in order to compute an approximation of the following expectation

    𝔼(f(Y(T)))=𝔼(f(∫_0^T e^μ (T-t)W_2(t) W(t)))

 for f=f(x) given as in (<ref>) with K=2. Computation of (<ref>) corresponds to derivative pricing, where the price of the underlying risky asset is described by (<ref>).
 
 The approximation to 𝔼(f(Y(T))) is defined by
 
    1/M∑_j=1^M f(𝒜_n,j^RM(X̃,W̃)),

 where M is a number of independent copies of  𝒜_n^RM(X̃,W̃). Due to the strong law of large numbers we get for all n∈ℕ

    1/M∑_j=1^M f(𝒜_n,j^RM(X̃,W̃))→𝔼(𝒜_n^RM(X̃,W̃)),   a.s.,

as M→ +∞. Moreover, since f:ℝ→ℝ is a Lipschitz function and X∈ F_L^1/2,2,q with q>2, the standard arguments and Theorem <ref> (i) yield the following estimate for averaged weak error, where  (X̃, W̃)∈ V_X(δ_1)×𝒲_s(δ_2) and δ_1,δ_2∈ [0,1],
 
    𝔼(f(ℐ(X, W)))-1/M∑_j=1^M f(𝒜_n,j^RM(X̃,W̃))_2   ≤ C_1·ℐ(X, W)-𝒜_n^RM(X̃,W̃)_2
       +2M^-1/2f(𝒜_n^RM(X̃,W̃))_2
       ≤ C_2(n^-1/2+δ_1+δ_2+M^-1/2).




 §.§ Noise

For the purpose of testing we analyze following disturbing functions

    p_1(t,x)    = 1, 
    
    	p_2(t,x)    = x, 
    
    	p_3(t,x)      = xt^2.


It is worth to mention that the noise function p_1 corresponds to the standard absolute deterministic noise and p_2, p_3 are related to the standard relative error. The latter can be connected with the computation precision. There is a trend now observed in computations, e.g. for deep learning, where the computations are conducted in lower precision in order to gain huge computation speedup. The novel GPU architectures (e.g. NVIDIA Volta) are designed with some dedicated accelerators for single or half precision operations. 

For each test the information about the analyzed precision level δ_1 and δ_2 will be given. All the tests were conducted with r=2. 



 §.§ Error criterion

For the problem (<ref>) we know the exact value of the solution, therefore we can have the following error estimate 


    ℐ(X_1,W)  - 𝒜_n^RM(X̃_1, W̃) _2 ≈( 1/M∑_j=1^M ( ℐ_j(X_1,W) - 𝒜_n,j^RM(X̃_1, W̃)  )^2)^1/2,

where M=2048 corresponds to the number of computed independent realizations under given precision levels. 
In case of the problems (<ref>)-(<ref>), the exact solution is not known, hence in order to analyze the algorithm error, we need to compare the obtained result with the result obtained on the same trajectories for denser mesh. In our tests, as the expected convergence ratio is of no less than 0.5, it is reasonable to have thousand times more points. That leads to following error estimation formula, used for (<ref>)-(<ref>)


    ℐ(X_m,W)  - 𝒜_n^RM(X̃_m, W̃) _2 ≈( 1/M∑_j=1^M ( 𝒜_L· n,j^RM(X_m, W) -𝒜_n,j^RM(X̃_m, W̃) )^2)^1/2,


where m∈{2,3,4} and L = 1000.  
For (<ref>) we use the following quantity
 
    |1/M∑_j=1^M f(𝒜_n,j^RM(X̃,W̃))-1/M∑_j=1^M f(𝒜_L· n,j^RM(X,W))|


  as the approximation of the weak error
  
    |𝔼(f(𝒜_n^RM(X̃, W̃)))-𝔼(f(Y(T)))|.

  



 §.§ Results

In  Figure <ref> we present the behavior of the error for the Riemann-Maruyama quadrature 𝒜^RM_n  for problem (<ref>). The numerical results are compared with the theoretical rate of convergence obtained for the algorithm, i.e. we present the effect of changing the precision levels δ_1 and δ_2. In Figures <ref>-<ref> we present behavior of the error for the Riemann-Maruyama quadrature 𝒜^RM_n for problems (<ref>)-(<ref>). (The  errors are measured accordingly to (<ref>).) From the Figure <ref> we see that if δ_1, δ_2 are on the level n^-1/2 then the  Riemann-Maruyama quadrature preserves the error 𝒪(n^-1/2), known from the case when the information is exact. The results confirm the necessity of tending with the precision parameters to zero in order to maintain the convergence rate for the Riemann-Maruyama quadrature.

Results for the weak approximation are given at Figure <ref>.



 §.§ GPU implementation

Below we present pseudo-code for the GPU implementation of the algorithm 𝒜_n^RM. This algorithm is designed for the case where we wish to compute multiple realizations of 𝒜_n^RM, returning the array of results. That algorithm, because of straightforward usage of parallel programming, enabled significant computational improvement of using graphics processing units. Moreover, additional speedup can be observed for using GPUs also for generating normally distributed numbers. Hence, it is suitable for e.g. derivative pricing. 







































In our experiments we compared the performance of the algorithm 𝒜_n^RM for both GPU and CPU implementations. For GPU implementation we used 32 blocks and 64 threads for problems (<ref>), (<ref>), (<ref>) and 512 threads for problem (<ref>). The CPU performance was tested using 8 threads. The used hardware was GPU – NVIDIA TESLA P100 (Maxwell), CPU – Intel Xeon E5-2680v4 (Broadwell). The speedup comparison for problem (<ref>) is presented in Figure <ref>. As we can see it is possible to have speedup of level 100x.












§ CONCLUSIONS
 
We investigated the problem of efficient approximation of Itô integrals under inexact information about the Wiener process and an integrand. We showed that for certain precisions (δ_1=δ_2≥ 0) the Riemann-Maryama quadrature rule is optimal. We also proposed  GPU implementation of the algorithm that is suitable for practical purposes.



Acknowledgments.


The authors were partially supported by the Faculty of Applied Mathematics AGH UST dean grant for PhD students and young researchers within subsidy of Ministry of Science and Higher Education, with grant numbers as follows: 15.11.420.038/18 (Andrzej Kałuża), 15.11.420.038/2 (Paweł M. Morkisz), and 15.11.420.038/1 (Paweł Przybyłowicz).

22
delmey
C. Dellacherie, P-A. Meyer, Probabilities and Potential, Hermann, 1978.

esikr
M. Eisenmann, R. Kruse, Two quadrature rules for stochastic Itô integrals with fractional Sobolev regularity, https://arxiv.org/abs/1712.08152.

Fatahalian
K. Fatahalian, J. Sugerman,  P. Hanrahan, Understanding the efficiency of GPU algorithms for matrix-matrix multiplication, Graphics Hardware (2005), 133–137.

alex
N. Whitehead, A. Fit-Florea, Precision & Performance: Floating Point and IEEE 754 Compliance for NVIDIA GPUs, NVIDIA, 2011.

legall
J-F. Le Gall, Brownian Motion, Martingales, and Stochastic Calculus, Springer, 2016.
 
hein1
S. Heinrich, Lower complexity bounds for parametric stochastic Itô integration, Preprint, https://www.uni-kl.de/AG-Heinrich/papers/lowpsint17.pdf,  2017. 

heindaun
S. Heinrich, T. Daun, Complexity of Banach space valued and parametric stochastic Itô integration, J. Complex. 40  (2017), 100–122.
 
hert
P. Hertling,  Nonlinear Lebesgue and Itô integration problems of high complexity, J. Complex. 17  (2001), 366–387.

 
JoC4
B. Kacewicz, M. Milanese, A. Vicino, Conditionally optimal algorithms and estimation of reduced order models,
J. Complexity 4 (1988), 73–85.


KaPl90
B. Kacewicz, L. Plaskota, On the minimal cost of approximating linear problems based on information with deterministic noise, Numer. Funct. Anal. and Optimiz. 11 (1990), 511-528.

KaPr16
B. Kacewicz, P. Przybyłowicz, On the optimal robust solution of IVPs with noisy information, Numer. Algor. 71 (2016),  505–518.

KARSHRI. Karatzas, S. E. Shreve, Brownian Motion and Stochastic Calculus, 2n edition, Springer - Verlag, New York, 1991.

KP P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag Berlin, Heidelber, third edition, 1999.

Kruger
J. Krüger, R. Westermann, Rüdiger, Linear Algebra Operators for GPU Implementation of Numerical Algorithms, ACM Trans. Graph. 22 (2003), 908–916.

Langdon W. B.  Langdon,  Brian Y. H. Lam, J. Petke, M.  Harman. Improving CUDA DNA Analysis Software with Genetic Programming. Proceedings of the 2015 on Genetic and Evolutionary Computation Conference - GECCO '15. (2015), 1063–1070.

NV_Volta NVIDIA, NVIDIA Tesla V100 GPU Architecture, available online http://images.nvidia.com/content/volta-architecture/pdf/volta-architecture-whitepaper.pdf.







 

MoPl16
 P. M. Morkisz, L. Plaskota, Approximation of piecewise Hölder functions from inexact information, J. Complex. 32 (2016), 122–136.
 
MoPr17 
P. M. Morkisz, P. Przybyłowicz, Optimal pointwise approximation of SDE's from inexact information, J. Comp. and Appl. Math. 324 (2017), 85–100.




Pla96
  L. Plaskota, Noisy Information and Computational Complexity, 
   Cambridge Univ. Press, Cambridge, 1996.

Pla14
  L. Plaskota, Noisy information: optimality, complexity, tractability, 
  in Monte Carlo and quasi-Monte Carlo Methods 2012, 
  J. Dick, F.Y. Kuo, G.W. Peters, I.H. Sloan (Eds.), Springer 2013,  173–209. 

PP1
P. Przybyłowicz, Linear information for approximation of the Itô integrals. 
Numer. Algorithms 52 (2009), 677–699. 


PP2 
P. Przybyłowicz, Adaptive Itô–Taylor algorithm can optimally approximate the Itô integrals 
of singular functions. J. Comp. Appl. Math. 235 (2010), 203–217.

PP3 
P. Przybyłowicz, Optimal sampling design for approximation of stochastic Itô integrals with 
application to the nonlinear Lebesgue integration. J. Comp. Appl. Math. 245 (2013), 10–29.

Riesinger
C. Riesinger, T. Neckel, F. Rupp, Solving random ordinary differential equations on GPU clusters using multiple levels of parallelism, SIAM J. Sci. Comput. 38(4) (2016), pp. C372–C402.

Riesinger2
C. Riesinger, T. Neckel, F. Rupp, A. P. Hinojosa, H.-J. Bungartz, GPU optimization of pseudo random number generators for random ordinary differential equations, Procedia Computer Sci 29 (2014), pp. 172-183.




Ryoo
S. Ryoo, C. I. Rodrigues, S. S. Baghsorkhi, S. S. Stone, D. B. Kirk, W. W. Hwu, Optimization principles and application performance evaluation of a multithreaded GPU using CUDA, Proceedings of the 13th ACM SIGPLAN PPoPP, (2008), 73-82.















TWW88
  J.F. Traub, G.W. Wasilkowski,  H. Woźniakowski,
  Information-Based Complexity, Academic Press, New York, 1988.

waswo
G. W. Wasilkowski, H.  Woźniakowski, On the complexity of stochastic integration, Math. Comp. 70 (2001), 685-698.


Wer96
  A.G. Werschulz, The complexity of definite elliptic problems with noisy data. J. Complex. 12 (1996), 440-473.
  
Wer97
  A.G. Werschulz, The complexity of indefinite elliptic problems with noisy data. J. Complex. 13 (1997), 457-479.


