Mean Field Game MURI Report Stanley Osher (UCLA) August, 2022 1 AFOSR MURI FA9550-18-1-0502

2

MURI’s Product 3

Mean fi eld game 4

Mean fi eld system A typical time-dependent MFG system (derived from (9)) has the form 8 > > > > > < > > > > > : t 2 + H (x, Dx (x, t), ⇢(x, t), t) = 0 ⇢t 2 ⇢ divx (⇢DpH (x, Dx (x, t), ⇢(x, t), t)) = 0 ⇢ 0, R Td ⇢(x, t)dx = 1, t 2 [0, T] ⇢(x, T) = ⇢T (x), (x, 0) = 0(x, ⇢(x, 0)), x 2 Td . (10) Given: H : Td ⇥ Rd ⇥ X ⇥ R ! R is a periodic in space and convex in momentum Hamiltonian, where X = R+ or X = L 1(Td ; R+) or X = R+ ⇥ L 1(Td ; R+). H = L ⇤ is a di↵usion parameter 0 : Td ⇥ X ! R, ⇢T : Td ! R given initial-terminal conditions Unknowns: , ⇢ : Td ⇥ [0, T] ! R has the form 5

Mean field information dynamics via reaction di↵usions Wuchen Li (USC), Wonjun Lee (UCLA), Stanley Osher (UCLA) The paper will be appeared in the Journal of Computational Physics. 1 6

Reaction di↵usion Consider a scalar nonlinear reaction di↵usion equation by @tu(t, x) = F(u(t, x)) + R(u(t, x)), where x 2 ⌦, ⌦ ⇢ Rd is a compact convex spatial domain, u 2 M(⌦) = {u 2 C1(⌦): u 0}, and is the Euclidean Laplacian operator. 2 7

Di↵usions and Entropy dissipation Consider a heat equation in Rd by @u(t, x) @t = r · (ru(t, x)) = u(t, x). Consider the negative Boltzmann-Shannon entropy by H(u) = Z Rd u(x)log u(x)dx. Along the heat equation, the following dissipation relation holds: d dtH(u(t, ·)) = Z Rd krx log u(t, x)k 2u(t, x)dx = I(u(t, ·)), where I(u) is named the Fisher information functional. There is a formulation behind this relation, namely I The heat flow is a gradient flow of entropy in optimal transport metric. 3 8

Gradient flows The Wasserstein gradient flow of an energy functional F(u) leads to @tu = G(u) 1 uF(u) =r · (ur uF(u)). I If F(u) = R F(x)u(x)dx, then @tu = r · (urF(x)). I If F(u) = R u(x) log u(x)dx, then @tu = r · (ur log u) = u. 4 9

Lyapunov methods induced metrics I Denote the smooth positive probability space by M = n u 2 C1(⌦): u > 0 o . Given F, G: R ! R satisfying F0(u) > 0 if u > 0, and G00(u) > 0. Denote V1(u) = F0(u) G00(u) , V2(u) = R(u) G0(u) . Denote the tangent space of M at u 2 M by TuM = n 2 C1(⌦) o . 5 10

Lyapunov methods induced metrics II Definition (Mean-field information metric) The inner product g(u): TuM ⇥ TuM ! R is given below. For any 1 , 2 2 TuM, define g(u)( 1, 2) = Z ⌦ 1 ⇣ r · (V1(u)r) + V2(u) ⌘ 1 2dx = Z ⌦ (r 1, r 2)V1(u)dx + Z ⌦ 1 2V2(u)dx, where i = r · (V1(u)r i) + V2(u) i, i = 1, 2. 6 11

Reaction-di↵usions in metric spaces Given an energy functional F : M ! R, the gradient flow of F in (M(⌦), g) satisfies @tu(t, x) = r · (V1(u)r uF(u))(t, x) V2(u) uF(u)(t, x). In particular, if F(u) = G(u) = R G(u)dx, then @tu =r · (V1(u)rG0(u)) V2(u)G0(u) =r · ( F0(u) G00(u) G00(u)ru) + R(u) G0(u) G0(u) =r · (F0(u)ru) + R(u) = F(u) + R(u). 7 12

Controlling reaction di↵usions Denote a given energy functional by F : M(⌦) ! R. Consider a variational problem by inf v1,v2,u Z 1 0 h Z ⌦ 1 2kv1k 2V1(u) + 1 2|v2| 2V2(u)dx F(u) i dt where the infimum is taken among all density functions u: [0, 1] ⇥ ⌦ ! R, vector fields v1 : [0, 1] ⇥ ⌦ ! Rd, and reaction rate functions v2 : [0, 1] ⇥ ⌦ ! R, such that @tu(t, x) + r · (V1(u(t, x))v1(t, x)) = v2(t, x)V2(u(t, x)), with fixed initial and terminal density functions u0 , u1 2 M(⌦). 8 13

Mean field information dynamics The critical point of variational problem satisfies v1(t, x) = r (t, x), v2(t, x) = (t, x), with 8 > > > > < > > > > : @tu(t, x) + r · (V1(u(t, x))r (t, x)) = V2(u(t, x)) (t, x) @t (t, x) + 1 2kr (t, x)k 2V 0 1 (u(t, x)) + 1 2| (t, x)| 2V 0 2 (u(t, x)) + uF(u)(t, x) = 0. 9 14

Example I Example (Wasserstein metric, heat flow and density Hamilton flows) Let G(u) = u log u u, F(u) = u, R(u) = 0, thus V1(u) = F0(u) G00(u) = u, V2(u) = R(u) G0(u) = 0. Then the mean-field information metric forms g(u)( 1, 2) = Z ⌦ (r 1(x), r 2(x))u(x)dx. Here the Hamiltonian flow in (M(⌦), g) satisfies 8 < : @tu + r · (ur ) = 0 @t + 1 2kr k 2 + uF(u) = 0. 10 15

Example II Example (Fisher-KPP metric, Fisher-KPP equation and density Hamiltonian flows) Consider F(u) = u, G(u) = u log u u, R(u) = u(1 u). Thus V1(u) = F0(u) G00(u) = u, V2(u) = R(u) G0(u) = u(u 1) log u . Here the Hamiltonian flow in (M(⌦), g) satisfies 8 > > < > > : @tu + r · (ur ) u(u 1) log u = 0 @t + 1 2kr k 2 + | | 2 (2u 1) log u + 1 u (log u)2 + uF(u) = 0. 11 16

G-prox Primal-dual methods Consider the following convex optimization problem. min z f(Az) + g(z), f and g are convex, We use G-Prox PDHG to solve the problem by iterating p (k+1) = argmax p L(z (k+1) , p) 1 2 kp p (k)k2 H z (k+1) = argmin z L(z, 2p (k+1) p (k)) + 1 2⌧ kz z (k)k2 L 2 . where p is a dual variable, L is a Lagrangian functional L(z, p) = g(z) + hAz, pi f ⇤(p) and f ⇤(p) = sup z hz, pi f(z) is the Legendre transform of f. G-Prox PDHG chooses a specific norm H which allows step sizes to be independent of an operator A. ⌧kAk2 < 1 (classical PDHG) ! ⌧ < 1 (G-Prox PDHG). 12 17

Algorithm Algorithm 2: G-prox PDHG for mean-field information variational problem Input: Initial density u0. Output: u, , m2 : [0, T] ⇥ ⌦ ! R, m1 : [0, T] ⇥ ⌦ ! Rd. For k 2 N For (t, x) 2 [0, 1] ⇥ ⌦ (k+1)(t, x) = (k)(t, x) + (AAT ) 1 ⇣ @ t u(k)(t, x) + r · m(k) 1 (t, x) m(k) 2 (t, x) ⌘ ; m(k+1) 1 (t, x) = V1(u (k)(t,x)) ⌧+V1(u (k)(t,x)) ⇣ m(k) 1 (t, x) + ⌧r 2 (k+1)(t, x) (k)(t, x) ⌘ ; m(k+1) 2 (t, x) = V2(u (k)(t,x)) ⌧+V2(u (k)(t,x)) ⇣ m(k) 2 (t, x) + ⌧ 2 (k+1)(t, x) (k)(t, x) ⌘ ; u(k+1)(t, x) Compute using Newton’s method u(k+1)(1, x) = (G⇤)0( 2 (k+1)(1, x) + (k)(1, x)). 13 18

Numerical Example We present two sets of numerical experiments using the Algorithm with di↵erent V1 , V2 functions. We used the discretization size 128 ⇥ 128 ⇥ 30 (x-axis ⇥ y-axis ⇥ time) In this example, we use V1(u) = u, V2(u) = u↵, ↵ > 0 We show how varying the exponent ↵ a↵ects the evolution of the densities. Below are the initial and terminal densities used in the simulations. u0(0, x) = 10 exp 60 (x 0.3)2 + (y 0.7)2 + 1 u1(0, x) = 20 exp 60 (x 0.7)2 + (y 0.3)2 + 1. 14 19

Figure: Example 1: Snapshots of the densities at di↵erent times t. The first, the second, the third and the fourth rows represent V2(u) = u, V2(u) = u 2 , V2(u) = u 3 , V2(u) = u 4 , respectively. 15 20

Figure: Example 1: Cross sections along a line x + y = 1 for di↵erent time t. Each color represents V2(u) = u, V2(u) = u 2 , V2(u) = u 4 . 16 21

A primal-dual approach for solving conservation laws with implicit in time approximations Siting Liu, Stanley J. Osher, Wuchen Li, and Chi-Wang Shu arXiv: 2207.07234 2022 22

Introduction We consider the following initial value problem (IVP) of scalar conservation law defined over the domain ⌦ ⇥ [0, T]: ( @tu(x, t) + @xf(u(x, t)) @x ( (x)@xu (x, t)) = 0, u(x, 0) = u0(x). Computational methods I High order accurate schemes in the spatial domain: discontinuous Galerkin methods (DG methods), ENO, WENO, etc. I Explicit high order schemes in the time domain: the forward Euler method and Runge-Kutta method, etc. They have restrictions on the length of the time step from Courant-Friedrichs-Lewy (CFL) conditions. I Implicit methods are suitable for sti↵ problems and can be used with larger time steps. The solvers typically require: extra computation and storage of the Jacobian matrix; a good initial guess etc. 2 23

A new framework for solving conservation laws We connect the computational methods with primal-dual optimization approach. I Cast IVP as a saddle point problem (a min-max problem). – Treat the continuity equation as a constraint. – Get a system of dual equations. I Apply implicit in time scheme for the min-max problem. – Accommodate many high order scheme in time and space. – Achieve implicit schemes for system of dual equations. I Use primal-dual method to find the saddle point. – First order optimization approach, no computation or storage of Jacobian matrix, no requirement on initial guess. – Optimization iterations independent of meshsize. – Geometric convergence for linear equations. No theoretic guarantee for other cases, but convergence in practice. – Highly parallelizable and easy-to-implement. No nonlinear inversions are required! 3 24

A saddle point formulation of the IVP We recast the IVP as the following optimization problem: min u 1u2U . U = n u : u(x, 0) = u0(x), @tu(x, t) + @xf(u(x, t)) @x ( (x)@xu (x, t)) = 0 for all (x, t) 2 ⌦ ⇥ [0, T] o . Here the indicator function defined as follows 1a2A = ( 0 if a 2 A 1 if a / 2 A. By introducing a Lagrange multiplier : ⌦ ⇥ [0, T] ! R, we can remove the constraint, which leads to a min-max problem: min u,u(·,0)=u0 max L(u, ), (1) L(u, ) = Z T 0 Z ⌦ (x, t) (@tu(x, t) + @xf(u(x, t)) @x ( (x)@xu (x, t))) dxdt. 4 25

The system of dual equations at saddle point We consider a variable coe cient type heat equation, then L(u, ) = Z T 0 Z ⌦ (ut @x ( (x)@xu)) dxdt = Z T 0 Z ⌦ u(x, t) (@t (x, t) + @x ( (x)@x )) dxdt + Z ⌦ u(x, T) (x, T) u(x, 0) (x, 0)dx. By taking the first-order optimality condition of L(u, ), we obtain the following system of equations: u ( ) runs forward (backward) in time: 8 > > > < > > > : ut @x ( (x)@xu (x, t)) = 0, u(x, 0) = u0(x), t + @x ( (x)@x (x, t)) = 0, (x, T) = 0. 5 26

Finite di↵erence scheme for variable coe cient type heat equation The discretized min-max problem (1) takes the following formulation: min u2{ul j } max 2{ l j } L(u, ), where L(u, ) = hthx X 0lNt 1 1jNx l j ul+1 j ul j ht Lap(ul+1)j ! . The discretized saddle point of the min-max problem is: 8 > > > > > > > < > > > > > > > : ul+1 j ul j ht Lap(ul+1)j = 0, 0  l  Nt 1, 1  j  Nx, u0 j = u0(xj), 1  j  Nx, l j l 1 j ht + Lap( l 1)j = 0, 1  l  Nt, 1  j  Nx, Nt j = 0, 1  j  Nx. Both equations are discretized with implicit backward Euler method. 6 27

A primal-dual algorithm We adapt the primal-dual hybird gradient (PDHG) algorithm for convex optimization problem, and integrate Gneral-prox PDHG and nonlinear PDHG. The n-th iteration is taken as follows: 8 > > > < > > > : un = argmin hu, AT ˜n 1i + 1 2⌧u ku un 1k2 L2 , n = argmax hAun, i 1 2⌧ k n 1k2 H , ˜n = 2 n n 1. Here A denotes the discrete approximation of the di↵erential operator: Au ⇡ @tu @x ( (x)@xu) , AT ⇡ t + @x ( (x)@x ) . Each update can be written explicitly, un = un 1 ⌧u ⇣ (@t + ↵@x) ˜n 1 + @x ⇣ @x ˜n 1 ⌘⌘ , n = n 1 ⌧ @tt + ↵2@xx ˆ2@xxxx 1 ((@t + ↵@x) un @x( @xun)) . The computation for n can be done by using Fast Fourier Transform (FFT). 7 28

Example I: variable coe cient type heat equation ( @tu = @x ( @xu) , (x, t) 2 [0, 1] ⇥ [0, 0.1], u(x, 0) = exp ( 64(x 0.5)2), for x 2 [0, 1], We record number of iterations needed to achieve certain errors ✏ for the primal-dual approach for di↵erent mesh in the following table: Nt ⇥ N 64 ⇥ 16 128 ⇥ 32 256 ⇥ 64 512 ⇥ 128 1028 ⇥ 256 ✏ = 10 6 90 91 92 94 102 ✏ = 10 10 144 146 147 149 163 I For an explicit scheme, we need to satisfy a restrictive CFL condition, i.e., ht (hx)2 < 1 2 maxx (x) . However, in the above examples we have ht (hx)2 = 128 5 , 256 5 , ..., 2048 5 . I We can see from above that the convergence of the algorithm for this linear equation is independent of the grid size. 8 29

Example II: quadratic conservation law Consider a quadratic conservation law of the following form: @tu + @x(↵u2 + u) = 0, for (x, t) 2 [0, 2] ⇥ [0, 1], u(x, 0) = ( 0.1, if 1  x  2, 0.25 else. For the numerical scheme, we use backward Euler and discontinuous Galerkin method. Figure: Here ↵ = 1, = 1, we have the tra c equation. Left: numerical solution of the tra c equation with discontinuous initial data, N ⇥ Nt = 256 ⇥ 32. Right: numerical solution at time t = 0, 0.5, 1. The shock propagation (left) and development of rarefaction wave (right). 9 30

Example III: Barotropic Compressible Navier–Stokes equation (BNS) Consider a BNS equation of the following form: 8 > < > : @t⇢ + @xm = 0, @tm + @x( m2 ⇢ ) + @xP(⇢) = @x(µ(⇢)@x m ⇢ ), Here, ⇢: density; m := ⇢v: flx; v: velocity. We set pressure P(⇢) = c1⇢2, µ(⇢) = c2 . For the numerical scheme, we use backward Euler and Lax–Friedrichs type of scheme. Figure: From left to right: ⇢, m, v. With positive uniform initial velocity start at t = 0, the density moves towards right. 10 31

Computational Mean-field Games on Manifolds Jiajia Yu (RPI) 1 Rongjie Lai (RPI) 1 Wuchen Li (U of SC) Stanley Osher (UCLA) June, 5th, 2022 1J. Yu and R. Lai’s work are supported in part by an NSF Career Award DMS-1752934 and NSF DMS-2134168 32

MFG on Manifolds I Conventional mean-field games deal with Nash Equilibrium on flat Euclidean domains. I But Euclidean domain is not enough for real-life applications, for example modeling evolution on the Earth or on data-sets or on the brain surface. I We generalize mean-field games to manifolds, and show that with operators on manifolds, the Nash Equilibrium satisfies 8 > < > : @t (x, t) + H(x, rM (x, t)) = F(x, ⇢(·, t)), @t⇢(x, t) rM · (⇢(x, t)@qH(x, rM (x, t)) = 0, (x, 1) = FT (x, ⇢(·, 1)), ⇢(·, 0) = ⇢0. Figure source: NASA(top), MNIST, Wikipedia(center), [Lai et al’13](bottom). 2 33

Variational Formulations I The PDE system is very complicated in its coordinate form 8 > > > > < > > > > : @t X (⇠ 1 , ⇠ 2 , t) + HX (X, g 1 ⇠ r⇠ X ) = F(X, ⇢(·, t)), @t⇢X (⇠ 1 , ⇠ 2 , t) 1 p det(g⇠ ) 2 X i=1 @ @⇠i 2 X j=1 p det(g⇠ )⇢X (g 1 ⇠ )ij@qj HX ! = 0, X (⇠ 1 , ⇠ 2 , 1) = FT (X(⇠ 1 , ⇠ 2), ⇢(·, 1)), ⇢X (⇠ 1 , ⇠ 2 , 0) = ⇢0(X(⇠ 1 , ⇠ 2)). I Comparing to the PDE formulation, when F(⇢) ⇢ (x) = F(x, ⇢), FT (⇢) ⇢ (x) = FT (x, ⇢), the equivalent variational formulation is easier to handle min ⇢,m Z 1 0 Z M ⇢(x, t)L ✓ x, m(x, t) ⇢(x, t) ◆ dMxdt + Z 1 0 F(⇢(·, t))dt + FT (⇢(·, 1)), s.t. @t⇢(x, t) + rM · m(x, t) = 0, ⇢(·, 0) = ⇢0. A local coordinate representation of a manifold. 3 34

Space Discretization I With the triangular mesh approximation of manifold and linear functions on the manifold, we discretize the problem to min P,M e Y(P, M) := Z 1 0 s X j=1 ATj P(Tj, t) kM(Tj, t)k2 2 dt + Z 1 0 e F(P(·, t))dt + e FT (P(·, 1)) s.t. A(P, M) := ✓ d dt P(Vi, t) + (r f M · M)(Vi, t) ◆ Vi2V,t2[0,1] = 0, P(·, 0) = P0. A kitten triangular mesh. A linear function on a mesh. The gradient on a triangle. 4 35

Algorithm I We applied projected gradient descent to solve the variational problem. 8 > > > < > > > : (P, M)(k+1/2) := ( b P, c M)(k) ⌘ (k)rP,M e Y( b P, c M)(k) , (P, M)(k+1) := proj{A(P,M)=0} (P, M)(k+1/2) , ( b P, c M) := ⇣ 1 + ! (k) ⌘ (P, M)(k+1) ! (k)(P, M)(k) . I The projection operator is exactly (P, M)(k+1) = (Id A⇤(AA⇤) 1A)(P, M)k+1/2 . I We pre-compute (AA⇤) 1 to save total computational cost. 5 36

Numerical Results: Avoiding Obstacles I Recall: the objective function is Z 1 0 Z M km(x, t)k2 g(x) 2⇢(x, t) dMxdt + Z 1 0 F(⇢(·, t))dt + FT (⇢(·, 1)) I The interaction cost F(⇢(·, t)) = R M ⇢(x, t)b(x, t)dMx encourages the density to avoid obstacles. I Our model and algorithm can handle spherical or genus-2 manifolds very well. US-map-based mesh. “8”-shape 6 37

Numerical Results: Concentrated Density I Recall: the objective function is Z 1 0 Z M km(x, t)k2 g(x) 2⇢(x, t) dMxdt + Z 1 0 F(⇢(·, t))dt + FT (⇢(·, 1)) I The interaction cost F(⇢(·, t)) = R M p ⇢(x, t) + ✏dMx encourages the density to concentrate. F(⇢(·, t)) = 0. F(⇢(·, t)) = R M p ⇢(x, t) + ✏dMx. 7 38

Numerical Results: Smooth Density I Recall: the objective function is Z 1 0 Z M km(x, t)k2 g(x) 2⇢(x, t) dMxdt + Z 1 0 F(⇢(·, t))dt + FT (⇢(·, 1)) I The interaction cost F(⇢(·, t)) = 1 2 R M krM⇢(x, t)k2 g(x) dMx encourages the density to be smooth. F(⇢(·, t)) = 0. F(⇢(·, t)) = 1 2 R M krM⇢(x, t)k2 g(x) dMx. 8 39