Particle algorithm for McKean SDE: a short review on numerical - PowerPoint PPT Presentation

Particle algorithm for McKean SDE: a short review on numerical analysis Mireille Bossy Sophia Antipolis CEMRACS 2017– Numerical methods for stochastic models July 17 - August 25, CIRM, Marseille

Let’s start with an example of simulation with particles Stochastic Lagrangian model for wind simulation, forced with a WRF simulation, in the Marseille navigation basin (60 km × 60 km × 2 km. 150 × 150 × 200 computation cells, with 96 particles in each cell ( 432 × 10 6 particles). Time step is 10 seconds, with screen-out each 600 seconds. Simulation realized in collaboration with Yann Amice (Metigate) - 2016 2

Outline � Stochastic particle methods for McKean SDEs Method and introduction to the numerical analysis, for smooth coefficients � When the interaction kernels are non smooth functions Motivated by applications Atmospheric boundary layer simulation Wind farm simulation � Numerical analysis and numerical experiments for Stochastic Lagrangian models (toy version) 3

Probabilistic Interpretation of PDEs Linear framework, forward parabolic PDE (in one dimension) A PDE of conservation form  2 σ 2 ∂ 2 ρ t ∂ρ t ∂x 2 − ∂  ∂t = 1 ∂x ( B ( t, x ) ρ t ) in R , (1)  ρ t =0 = ρ 0 . B ( t, x ) , bounded measurable (1) is a Fokker-Planck Equation. � Set ( µ, f ) = R f ( x ) µ ( dx ) . ρ is a weak solution of (1), if for any test function f ∈ C 2 b ( R ) . � t ( ρ s , 1 2 σ 2 f ′′ + B ( s, · ) f ′ ) ds. ( ρ t , f ) = ( ρ 0 , f ) + (2) 0 When ρ 0 is a probability measure, (1) (or (2)) describes the dynamics of the time-marginals of a stochastic process. 4

Associated martingale problem Let ( X t ) t ≥ 0 the canonical process on ( C ([0 , + ∞ ) , R ) , B ( C ([0 , + ∞ ) , R ))) . For a probability P on ( C ([0 , + ∞ ) , R ) , B ( C ([0 , + ∞ ) , R ))) , we set P t = P ◦ X − 1 . t Definition - Martingale Problem P on ( C ([0 , + ∞ ) , R ) , B ( C ([0 , + ∞ ) , R ))) is a solution of the martingale problem 2 σ 2 ∂ 2 associated to L = 1 ∂x 2 + B ( t, x ) ∂ ∂x if i) P 0 = ρ 0 , ii) ∀ f ∈ C 2 b ( R ) , the process � t � 1 � 2 σ 2 f ′′ ( X s ) + B ( s, X s ) f ′ ( X s ) f ( X t ) − f ( X 0 ) − ds is a P -martingale. 0 5

Existence and uniqueness for X solution of � t X t = X 0 + B ( s, X s ) ds + σW t . 0 For any test function f � t � 1 � 2 σ 2 f ′′ ( X s ) + B ( s, X s ) f ′ ( X s ) E f ( X t ) = E f ( X 0 ) + ds. E (3) 0 By setting ρ t = P ◦ X − 1 the t -time-marginal of P , t � E f ( X t ) = R f ( x ) ρ t ( dx ) = ( ρ t , f ) and (3) rewrites (2) : � t ( ρ s , 1 2 σ 2 f ′′ + B ( s, · ) f ′ ) ds. ( ρ t , f ) = ( ρ 0 , f ) + 0 6

Numerical approximation, particle method Main idea : approximate the measure ρ t by the empirical measure of a set of N independent trials of the r.v. X t . Introducing N particles of independent dynamics � t X i t = X i B ( s, X i s ) ds + σW i 0 + t , i = 1 , . . . , N, 0 where (( X i 0 ) i =1 ,...,N ) are i.i.d with law ρ 0 , independent of the family of independent Brownian motions ( ( W i t ) , i = 1 , . . . , N ) . N t = 1 � N U δ X i t . N i =1 N Do we have the convergence of the empirical measure U t toward ρ t ? 7

(Basic) numerical analysis • The parabolic problem has a smooth solution ∀ t > 0 , ρ t ( dx ) = ρ t ( x ) dx. Introduce a smoothing empirical measure by a cut–off function N � � ( x ) = 1 � N,ε N φ ε ( x − X i ( x ) = φ ε ∗ U t ) U t t N i =1 � where φ ε ( x ) = 1 ε φ ( x ε ) and φ : R → R smooth and such that R φ ( x ) dx = 1 . At a given time t , we want the rate of convergence of the error : � � N,ε � � � ρ t ( x ) − U sup E ( x ) � . t x ∈ R 8

Set ρ ε t ( x ) = ( ρ t ∗ φ ε ) ( x ) . � � | x | θ | φ ( x ) | dx < ∞ ) x q φ ( x ) dx = 0 , ∀ q ≤ θ − 1 , If φ is a θ -order cut–off ( R R and if ρ t is smooth enough � � ∂ θ ρ t � � � ρ ε t ( x ) − ρ t ( x ) � L ∞ ( R ) ≤ Cε θ � � � � ∂x θ � � L ∞ ( R ) and � � � � � ≤ Cε θ + sup N,ε N,ε � ρ ε � � � � sup � ρ ( t, x ) − U ( x ) t ( x ) − U ( x ) E E t t � x ∈ R x ∈ R with � ρ ε t ( x ) = φ ε ( x − y ) ρ t ( y ) dy = E φ ε ( x − X t ) R N ( x ) = 1 � N,ε φ ε ( x − X i t ) . U t N i =1 The convergence is ensured by the strong law of large numbers. 9

More precisely, � �� N � � � � � 1 � N,ε � ( x ) − ρ ε � � φ ε ( x − X i � t ) − E φ ε ( x − X t ) E � U t ( x ) � = E t � � N � � i =1 � 1 E ( φ ε ( x − X t ) − E φ ε ( x − X t )) 2 ≤ √ N � 1 E ( φ ε ( x − X t )) 2 ≤ √ N 1 1 1 � ≤ √ √ ε � φ � 2 E φ ε ( x − X t ) L ∞ ( R ) N 1 1 1 � E ρ t ( Y ε − x ) = √ √ ε � φ � 2 L ∞ ( R ) N 1 1 1 ≤ √ � φ � L ∞ ( R ) � ρ t � 2 2 L ∞ ( R ) . εN 10

Adding a time discretisation of X i i i i k ∆ t ) + σ ( W i ( k +1)∆ t − W i X ( k +1)∆ t = X k ∆ t + ∆ tB ( k ∆ t, X k ∆ t ) . N � � � � 1 � N,ε N,ε, ∆ t i � � � φ ε ( x − X i � � k ∆ t ( x ) − U ≤ k ∆ t ) − φ ε ( x − X E � U ( x ) E k ∆ t ) k ∆ t � � N i =1 � � � � 1 � φ ′ � X 1 � E � � ≤ sup ε ( x ) k ∆ t − X t ∆ t � x ∈ R C ∆ t ε 2 (could by reduced to C ∆ t ≤ ε ) . Lemma � � � � N,ε, ∆ t N,ε, ∆ t � � � � sup � u ( t, x ) − U ( x ) � + E � u ( t, x ) − U ( x ) E t t � L 1 ( R ) x ∈ R C 2 + C 3 ∆ t ≤ C 1 ε θ + √ ε . εN 11

McKean Vlasov Fokker-Planck prototype  d d ∂ 2 � � ∂ ∂ρ t  ∂t = 1  ∂x i ∂x j ( a ij [ x, ρ t ] ρ t ) − ∂x i ( b i [ x, ρ t ] ρ t ) 2 (4) i,j =1 i =1   ρ 0 a given probability measure on R d . where, for a probability measure u � R d b ( x, y ) u ( dy ) ; for b ( x, y ) R d -valued, b [ x, u ] = a [ x, u ] = σ [ x, u ] t σ [ x, u ] , � σ [ x, u ] = R d σ ( x, y ) u ( dy ) ; σ ( x, y ) ( d × k ) matrix-valued. Distribution equation : for any test function f ∈ C 2 b ( R d ) , � � � t d d ∂ 2 f ρ s , 1 � � b i [ x, ρ s ] ∂f ( ρ t , f ) = ( ρ 0 , f ) + a ij [ x, ρ s ] ∂x i ∂x j + ds (5) 2 ∂x i 0 i,j =1 i =1 12

Probabilistic point of view Theorem [Méléard 95] for b and σ Lipschitz on R 2 d , for ρ 0 admitting a second order moment, there is existence and trajectorial uniqueness (given X 0 ∈ L 2 (Ω) and ( W t ) ) and in law of the solution to the SDE � t � t X t = X 0 + σ [ X s , ρ s ] dW s + b [ X s , ρ s ] ds, 0 0 where ρ s is the law of X s . (Proof from Sznitman 89, with σ = 1 , and bounded b ). P 2 = { probability P on C ([0 , T ] , R d ) such that E P (sup t ∈ [0 ,T ] | X 2 t | ) < ∞} , endowed with the Wasserstein metric. Φ : P 2 ∋ m �→ Law ( X m t , t ∈ [0 , T ]) ∈ P 2 � t � t X m σ [ X m b [ X m = X 0 + s , m s ] dW s + s , m s ] ds. t 0 0 admits a fix point in P 2 . 13

Associated particle system N � Idea : replace ρ t = Law ( X t ) by the empirical measure µ N = 1 δ X i,N . N i =1 � t � t X i,N = X i σ [ X i,N , µ N s ] dW i b [ X i,N , µ N 0 + s + s ] ds, i = 1 , . . . , N. (6) t s s 0 0 where (( X i 0 ) i =1 ,...,N ) are i.i.d with law ρ 0 , independent of the family of independent Brownian motions ( ( W i t ) , i = 1 , . . . , N ) . (6) is well-posed when b ( x, y ) σ ( x, y ) are Lipschitz. What convergence of µ N t toward ρ t ? 14

Propagation of chaos (Sznitman 89) Let E be a separable metric space and ν a probability measure on E . A sequence of symmetric probabilities ν N on E N is ν -chaotic if for any φ 1 , . . . , φ k ∈ C b ( E ; R ) , k ≥ 1 , k � � � ν N , φ 1 ⊗ . . . ⊗ φ k ⊗ 1 . . . ⊗ 1 � ν, φ l � . lim = N →∞ l =1 Application : Let ( X 1 ,N , . . . , X N,N ) , the set of N particles satisfying (6), valued in ( C ([0 , T ]; R d )) N , with law P N . Theorem (Méléard 95) ) is ρ 0 -chaotic, the chaos propagates and P N is P -chaotic If the initial of ( X 1 ,N , . . . , X N,N 0 0 where P is the unique law of the McKean SDE. 15

The tow following assertions are equivalent : (Tanaka 82) • P N is P -chaotic • µ N converges weakly toward the (deterministic) measure P Theorem (Méléard 95) � � � � | X i,N | 2 | X i t | 2 sup sup + sup sup < ∞ E E t N 0 ≤ t ≤ T N 0 ≤ t ≤ T � � | X i,N − X i t | 2 sup N E sup < ∞ , t N 0 ≤ t ≤ T where � t � t X i t = X i σ [ X i s , ρ s ] dW i b [ X i 0 + s + s , ρ s ] ds 0 0 with ρ t = P ◦ X − 1 . t 16

Numerical approximation � t � t X t = X 0 + σ [ X s , ρ s ] dW s + b [ X s , ρ s ] ds, with ρ s = Law ( X s ) . 0 0 Spatial discretisation : � t � t X i,N X i σ [ X i,N , µ N s ] dW i b [ X i,N , µ N s ] ds, i ≤ N. = 0 + s + s s t 0 ) Time-Discretisation : ∆ t > 0 , such that T = K ∆ t for K ∈ N ,  N � k ∆ t + ( 1 i,N j,N  X i,N ¯ ( k +1)∆ t = ¯ X i,N k ∆ t ))( W i ( k +1)∆ t − W i  σ ( X k ∆ t , X k ∆ t )  N    j =1  N + ∆ t ( 1 � b ( ¯ k ∆ t , ¯ X i,N X j,N k ∆ t )) ,    N   j =1   ¯ X i,N t =0 = X i 0 . N ( x ) = 1 � N,ε, ∆ t φ ε ( x − ¯ X i,N U k ∆ t ) , pour un cut-off donné φ . k ∆ t N i =1 17