Cubature Methods and Applications Dan Crisan Imperial College - PowerPoint PPT Presentation

Cubature Methods and Applications Dan Crisan Imperial College London CEMRACS 2017 Numerical methods for stochastic models: control, uncertainty quantification, mean-field July 17 - August 25 CIRM, Marseille Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 1 / 40

Talk Synopsis Feynman-Kac representations - a common platform Theoretical Analysis of solution of PDEs Approximation of Feynman-Kac representations Wiener measure approximation Computational effort Example 1: Semilinear PDEs (joint work with K. Manolarakis) The Feynman-Kac representation Functional discretization Theoretical results Numerical Implementation Example 2: linear parabolic SPDEs (joint work with S. Ortiz-Latorre) The Feynman-Kac representation Functional discretization Theoretical results Numerical Implementation Example 3: McKean-Vlasov PDEs (joint work with E. McMurray) The Feynman-Kac representation Functional discretization Theoretical results Numerical Implementation Final remarks Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 2 / 40

Feynman-Kac representations A common platform Common feature of many PDEs: their solutions can be represented as integrals of certain nonlinear functionals with respect to the Wiener measure. Feynman-Kac formula � u ( t , x ) = E [Λ t , x ( W )] = Λ t , x ( ω ) dP W ( ω ) ω ∈ C ([ 0 , ∞ ) , R d ) Microscopic level Macroscopic level Timeline � = 1 2 ∆ ut ∂ t ut Heat equation Feynman 1948 Kac 1949 u 0 = Φ Brownian motion Zakai equation dut = Lut + hut dYt Duncan, Mortensen, Zakai 1970 W = { Wt , t ≥ 0 } � d ∂ t ut = i , j = 1 aij ( ut ) ∂ i ∂ j ut McKean-Vlasov PDEs G¨ artner 1988 + � d i = 1 bi ( ut ) ∂ i ut + c ( ut ) ut � ∂ t ut = Lut + f � t , x , ut , ∇ ut � Semilinear PDEs Pardoux & Peng 1990, 1992 u 0 = Φ Fully Nonlinear PDEs F ( t , x , ut , ∇ ut , ∆ ut ) = 0 Soner, Touzi & Victoir 2007 3 − d incompressible � ∂ t ut + ( ut · ∇ ) ut − ν ∆ ut + ∇ p = 0 Constantin & Iyer 2008 Navier − Stokes equation ∇ · ut = 0 dut = Lut + ut (¯ K-S equation h )( dYt − ut ( h ) dt ) Crisan & Xiong 2009 ∂ t ut + ut ∂ x ut − ν∂ 2 viscous Burgers equation x ut = 0 Novikov & Iyer 2010 Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 3 / 40

Feynman-Kac representations Theoretical Analysis of Feynman-Kac formulae Smoothness of u t ≡ Smoothness of Λ t , x in Malliavin sense Ex: Heat equation � = 1 ∂ t u t 2 ∆ u t � 1 − ( y − x ) 2 u t ( x ) = E [Φ( x + W t )] = Φ( y ) √ e dy 2 t u 0 = Φ 2 π t Via an integration by parts formula one can prove that � u t ( x ) ′ = 1 ⇒ | u t ( x ) ′ | ≤ E [ | W t | ] 2 � φ � ∞ t E [ φ ( x + W t ) W t ] � φ � ∞ = √ . t π t Remarks: • Fundamental progress (though not complete) for F-K formulae for linear PDEs - notably through the Kusuoka-Stroock programme : Kusuoka & Stroock [1985,1987, 2003], further developed in DC & Ghazali [2007], DC, Manolarakis, Nee [2013], DC & Ottobre [2016]. • Some progress for F-K formulae for non-linear PDEs Semilinear equations : DC-Delarue [2012] McKean-Vlasov equations: DC & McMurray [2017] Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 4 / 40

Feynman-Kac representations Approximation of Feynman-Kac representation � u ( t , x ) = E [Λ t , x ( W )] = Λ t , x ( ω ) dP W ( ω ) ω ∈ C ([ 0 , ∞ ) , R d ) A three-step scheme: � n i = 1 δ ω i - ˜ W = 1 replace P W with P ˜ W approximates the signature of W n approximate Λ t , x with an explicit/simple version ˜ Λ t , x control the computational effort (use the TBBA) n u ( t , x ) ≃ 1 � ˜ Λ t , x ( ω i ) n i = 1 Full DNA Typical Paths Representative Paths Truncated DNA Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 5 / 40

Feynman-Kac representations Wiener measure approximation Chen’s iterated integrals expansion - signature of a path ⊗ i and T ( m ) � i = 0 ( R d ) ⊗ i be the tensor algebra = � ∞ = � m � R d � i = 0 ( R d ) R d � Let T of all non-commutative polynomials over R d and, respectively the tensor algebra of all non-commutative polynomials of degree less than m + 1. For a path ω ∈ C bv ([ 0 , ∞ ) , R d ) define its signature S t ( ω ) and, respectively, its truncated signature S m t ( ω ) to be the corresponding Chen’s iterated integrals expansion: ∞ � � R d � C bv ([ 0 , ∞ ) , R d ) → T � S : S t ( ω ) = d ω t 1 ⊗ ... ⊗ d ω t k 0 < t 1 ... t k < t k = 0 m � C bv ([ 0 , ∞ ) , R d ) → T ( m ) � R d � � S m S m d ω t 1 ⊗ ... ⊗ d ω t k . : t ( ω ) = 0 < t 1 ... t k < t k = 0 The (random) signature and, respectively, the truncated signature of the Brownian motion are ∞ m � � � � S m S t ( W ) = dW t 1 ⊗ ... ⊗ dW t k , t ( W ) = dW t 1 ⊗ ... ⊗ dW t k . 0 < t 1 ... t k < t 0 < t 1 ... t k < t k = 0 k = 0 Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 6 / 40

Feynman-Kac representations Wiener measure approximation E [ S t ( W )] uniquely identifies the Wiener measure P W . � � If ˜ t ( ˜ W is another process such that E [ S m S m t ( W )] = E W ) , then E [Λ t , x ( W )] ≃ E [Λ t , x ( ˜ W )] high order approximation of u ( t , x ) . See DC and Ghazali [2007] for conditions. Several choices for ˜ W : Kusuoka [2001,2004], Kusuoka and Ninomiya [2004], Lyons and Victoir. [2004], Ninomiya and Victoir [2004], etc. Theorem (Lyons & Victoir (2004)) For any t > 0 , there exists paths ω 1 , . . . , ω N ∈ C 0 0 , bv ([ 0 , t ]; R d ) and λ 1 , λ 2 , . . . , λ N ( � N i = 1 λ i = 1 ), such that if P ( ˜ W = ω i ) = λ i then � � t ( ˜ E [ S m S m t ( W )] = E W ) . W = � N If the above is true, we call L ˜ i = 1 λ i δ ω i the cubature measure and denote it by Q m t . Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 7 / 40

Feynman-Kac representations Wiener measure approximation For example, if we want to approximate E [ α ( X t )] , where X is the the solution of the following SDE d � V i ( X t ) ◦ dW i dX t = V 0 ( X t ) dt + t i = 1 Then X can be expressed as X = Ψ t ( W ) giving a representation of the form E [Λ t ( W )] . Choose X j to be the solution of the following ODE d dX j t = V 0 ( X j � V i ( X j t ) d ω j , i t ) dt + t i = 1 In this case: N E [Λ t ( ˜ � λ i g ( X j W )] = E Q m [ α ( X t )] = t ) i = 1 and m − 1 | E [Λ t ( W )] − E [Λ t ( ˜ 2 . W )] | ≤ C δ Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 8 / 40

Feynman-Kac representations Wiener measure approximation Cubature of order 3: For d=1, we can use 2 paths with equal weights λ j = 1 2 defined as ω j t = tz i , � d ( d + 2 ) � where z i ∈ {− 1 , 1 } . For d ≥ 2 , we can use + 1 linear paths with 6 equal weights λ j . Cubature of order 5: For d=1, we can use 3 paths, ω , − ω and 0 with ω is defined as √ √  � � 3 0 , 1 4 − t ∈ � � 22 t  2 3  √ √ √  √ � 1  � � � � � 3 t − 1 � 3 , 2 � ω ( t ) = 4 − + − 1 + t ∈ 22 3 22 6 3 3 √ √ √ √ � 2  � � � � �  3 3 t − 2 � � 2 + + 4 − t ∈ 3 , 1 22 22   6 2 3 Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 9 / 40

Feynman-Kac representations Computational effort reduction • An additional procedure is used to control the computational effort. • The measure Q m is replaced by a measure ˜ Q m , N with support of size N by using a tree based branching algorithm (TBBA). • The TBBA was introduced in DC & Lyons (2002) as the optimal stratified sampling procedure in the context of the filtering problem. The method has a wider applicability: it is applicable to any method that uses branching trees. • By merging the TBBA with the cubature method we keep the number of particles on the support of the intermediate measures constant . Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 10 / 40

Feynman-Kac representations Computational effort reduction • Assume that we constructed Q m = � M j = 1 λ j δ γ j with M paths and we want to reduce the number to at most N paths. • We replace Q m with a random measure ˜ Q m , N such that supp (˜ Q m , N ) ⊆ supp ( Q m ) and that the size on its support is at most N . For an arbitrary path γ ∈ supp ( Q m ) , we will have � ⌊ N Q m ( γ ) ⌋ with probability 1 − { N Q m ( γ ) } ˆ Q m N k ( γ ) = . (1) ⌊ N Q m ( γ ) ⌋ + 1 with probability { N Q m ( γ ) } N Q m , N is constructed so that it is a (random) probability measure, • In addition, ˜ i.e., � Q m ˆ k ( γ ) = 1 . (2) γ ∈ supp Q m • The mass allocated to each path γ ∈ supp ( Q m ) is either 0 or an integer Q m , N has size at most N . multiple of 1 / N ⇒ the support of any realization of ˜ Q m , N puts mass 1 / N on N • The maximum number of paths is achieved when ˜ of the M paths. This is the basis of the control of the computational effort. Dan Crisan (Imperial College London) Particle approximations to PDEs ans SPDEs 18 July 2017 11 / 40