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 4
inf ⇢,v J MFG (v, ⇢) = Z T 0 Z I R d L(x, v)⇢(x, t)dxdt + Z T 0 F(⇢(·, t))dt + G(⇢(·, T)), s.t. @ t ⇢(x, t) + r · (⇢(x, t)v(x, t)) = 0, ⇢(x, 0) = ⇢ 0 (x) (1) Our Contribution: solve these types of problems in high dimensions Our Approach I use Lagrangian coordinates to discretize the continuity equation I parameterize potential function (where v = r ) by a neural network 2 6
T) + c F (x, T) + G(z(x, T)) subject to @ t 0 B B B B @ z(x, t) l(x, t) c L (x, t) c F (x, t) c HJB (x, t) 1 C C C C A = 0 B B B B B @ r (z(x, t), t; ✓) (z(x, t), t; ✓) L ⇣ z(x, t), r (z(x, t), t; ✓) ⌘ F(z(x, t), t) c HJB (z(x, t), t) 1 C C C C C A , (2) where z(x, 0) = x, and the rest initialized with zero. I completely grid free I Laplacian of used to estimate determinant of induced transformation I (z(x, t), t; ✓) parameterized by neural network with weights ✓ 3 7
w>N(s, ✓ N )+ 1 2 s> ⇣ A + A> ⌘ s + c>s + b, ✓ = (w, ✓ N , vec(A), c, b), (3) where s = (x, t) 2 Rd+1 is the input feature and N(s, ✓ N ): Rd+1 ! Rm is a ResNet: u 0 = (K 0 s + b 0 ) u 1 = u 0 + h (K 1 u 0 + b 1 ) . . . . . . u M = u M 1 + h (K M u M 1 + b M ), (4) 4 8
computed as (s, ✓) = tr ⇣ E>(r2 s (N(s, ✓ N )w) + (A + A>))E ⌘ , where the columns of E 2 IRd+1⇥d are given by the ﬁrst d standard basis vectors in IRd+1. 5 9
computed as (s, ✓) = tr ⇣ E>(r2 s (N(s, ✓ N )w) + (A + A>))E ⌘ , where the columns of E 2 IRd+1⇥d are given by the ﬁrst d standard basis vectors in IRd+1. The trace of the ﬁrst layer of N(s, ✓ N ) can be computed as t 0 = ( 00(K 0 s + b 0 ) z 1 )>((K 0 E) (K 0 E))1, and we continue with the remaining rows in reverse order to obtain (N(s, ✓ N )w) = t 0 + h M X i=1 t i , t i = ( 00(K i u i 1 + b i ) z i+1 )>((K i J i 1 ) (K i J i 1 ))1, where J i 1 = r s u> i 1 2 IRm⇥d 5 10
⇢ 1 , target density d = 2, pull back d = 2, push fwd d = 2, characteristics d = 10, pull back d = 10, push fwd d = 10, characteristics d = 50, pull back d = 50, push fwd d = 50, characteristics d = 100, pull back d = 100, push fwd d = 100, characteristics qualitatively similar results in all dimensions 6 11
, target density pull back Lagrangian, ML push forward characteristics Eulerian, ﬁnite volume ⇢ 0 , initial density ⇢ 1 , target density pull back Lagrangian, ML push forward characteristics Eulerian, FV Haber E, Horesh, R, 2015 8 13
equilibrium solution can be found by minimizing an overall “energy” (e.g. multiply the value function for a single agent by ⇢): A(⇢, v) = min ⇢,v Z T 0 Z ⌦ ⇢(x, t)L(v) + F(x, ⇢) dx dt + Z ⌦ (x)⇢(x, T) dx where x = x(t), v = v(t) above, and where the optimization has the constraint, @t⇢ ⌫ ⇢ + div(⇢v) = 0. 3 16
objective, perform integration by parts, and push the minimization of v inside to obtain, = max ' min ⇢ Z T 0 Z ⌦ F(x, ⇢(x, t)) dx dt + Z ⌦ ⇣ (x) '(x, T) ⌘ ⇢(x, T) dx + Z ⌦ '(x, 0)⇢(x, 0) dx + Z T 0 Z ⌦ ✓ @t'(x, t) + '(x, t) L⇤( r'(x, t)) ◆ ⇢(x, t) dx dt which means we end up with a sampling problem. Then the idea is to turn ⇢ and ' into neural networks and train as in GANs (Generative Adversarial Networks). 4 17
D E x⇠real [D(x)] E z⇠N [D(G(z))] s.t.krDk 1 = min G max D E x⇠real [D(x)] E y⇠G(z),z⇠N [D(y)] s.t.krDk 1 I Training has a discriminator and a generator. The generator produces samples (analogous to our ⇢), and the discriminator evaluates the quality of those samples (analogous to our '). 5 18
turned into the min-max problem: max ' min ⇢ Z ⌦ ⇣ (x) '(x, T) ⌘ ⇢(x, T) dx + Z ⌦ '(x, 0)⇢(x, 0) dx + Z T 0 Z ⌦ F(x, ⇢(x, t)) dx dt + Z T 0 Z ⌦ ✓ @t'(x, t) + '(x, t) L⇤( r'(x, t)) ◆ ⇢(x, t) dx dt And supposing a certain form for F, then this can be expressed in expectation form, max ' min ⇢ E y⇠⇢(x,T ) [ (y) '(y, T)] + E y⇠⇢(x,0) ['(y, 0)] + Z T 0 E y⇠⇢(x,t) [F(y)] dt + Z T 0 E y⇠⇢(x,t) [@t'(y, t) + '(y, t) L ⇤( r'(y, t))] dt 6 19
is ⌫ = 0 (deterministic), and right is ⌫ = 0.4 (stochastic) I 100 dimensional obstacle problem with an obstacle in the ﬁrst two dimensions (“cylindrical” obstacle). 8 21
I 100 dimensional congestion problem with congestion in the ﬁrst two dimensions. ⌫ = 0 (left), and ⌫ = 0.1 (right) I Congestion term G(x, ⇢(x, t)) = R ⌦ R ⌦ 1 kx yk2+1 ⇢(x, T)⇢(y, T) dx dy, 9 22 G(T, ⇢)
to control-a ne systems. I Established existence results for the dynamical formulation for general Lagrangian costs and control constraints. I Numerically solved the optimal transport problem for the case of quadratic cost using Primal Dual Hybrid Gradient Method. Considered the homogeneous of degree one cost for the driftless systems. I Numerically solved the case of homogeneous of degree one cost by eliminating the time variable and hence reducing dimensionality of the problem. 2 24
robot modeled according to following di↵erential equation: 8 > < > : ˙ ✓ = u1 ˙ x1 = u2 cos ✓ ˙ x2 = u2 sin ✓, u1 : Turning Rate, u2 : Forward/Backward Velocity. Technical Challenge : Number of controls < Dimension of state space. 3 25
0 Z Rd L(t, x, u(t, x))⇢(t, x)dxdt subject to continuity equation constraint @t⇢ + r · (g0⇢) + n X i=1 r · (uigi⇢) = 0, initial and terminal constraints on density ⇢(0, ·) = ⇢0, ⇢(T, ·) = ⇢T . g0 : drift vector-ﬁeld, gi : control vector-ﬁelds. Goal: Find optimal control in feedback form u(t, x) such that solution ⇢(t, x) satisﬁes initial and terminal constraints and minimizes expected value of control Z T 0 Z Rd L(t, x, u(t, x))⇢(t, x)dxdt 4 26
ODE system ˙ x = g0(x) + n X i=1 vi(t)gi(x) and continuity, coercivity and convexity of the running cost L(t, x, u). I Establish existence of solutions for the optimal transport problem. I Extend previous theoretical results that only considered linear control systems and quadratic running cost L(t, x, u) = |u|2. 5 27
cost Wp (⇢0, ⇢T ), where p = 1, 2. Wp(⇢0, ⇢T )p = inf u,⇢ n Z T 0 Z ⌦ 1 p n X i=1 ⇢(t, x)|ui(t, x)|pdxdt s.t. CE holds. o When p = 2, we formalize an unconstrained min-max problem inf m,⇢ sup L(m, ⇢, ), L(m, ⇢, ) = Z Rd ⇢T (x) (T, x) ⇢0(x) (0, x)dx + Z T 0 Z Rd n X i=1 |mi(t, x)|2 2⇢ ⇢ ( t + r · g0(t, x)) r · n X i=1 mi(t, x)gi(x) ! dxdt. We apply G-prox Primal-Dual Hybrid Gradient Algorithm (PDHG) to solve the saddle point problem on the discretized domain. The rate of convergence independent of grid size, which is important for large scale problems. 6 28
model for two scenarios: initial density (green) splits into two parts (yellow) with di↵erent orientations. Case 1 (left) and Case 2 (right). The red arrows indicate the ﬁrst coordinate ✓ . Recall the control rule: 8 > < > : ˙ ✓ = u1 ˙ x1 = u2 cos ✓ ˙ x2 = u2 sin ✓, u1 : Turning Rate, u2 : Forward/Backward Velocity. 7 29
qualitatively di↵erent controls. For case 1, the densities steer and then split to their ﬁnal positions. For case 2, the density reach x1 -direction position ﬁrst, and then split to move horizontally towards ﬁnal positions. 8 30
cost with driftless system, we formalize the optimization problem to be independent of time-variable which reduces the dimensionality of the problem signiﬁcantly. 9 31
system with mixed coupling (nonlocal and local interactions). I Demonstrate that our approach is compatible with existing methods to solve single-agent control problems. I Combine ideas from convex optimization, kernel methods in machine learning and variational MFG. I Solve a class of non-potential MFG by casting it as primal-dual paris of monotone inclusions. 2 33
i ixiyi (Maximal spread) K2(x, y) = Y i µi exp ✓ |xi yi |2 2 2 i ◆ (Gaussian Repulsion) Nonpotential MFG: , + (x) = e x2 2 2 x<0 + e x2 2 2 + x 0 K3(x, y) = µ Y i i, , i,+ (yi xi) K4(x, y) = K3(Q 1x, Q 1y), Q : coordinates transformation matrix Feature expansions: K(x, y) = X i,j ki,j⇣i(x)⇣j(x) ⇣j(x): feature functions that capture the interactions between agents, e.g. Fourier series, polynomials etc. This is a computational e cient way of modeling the nonlocal interactions. 4 35
i exp ✓ |xi yi |2 2 2 i ◆ Smaller i yields stronger repulsion in xi direction. Obstacles modeled through soft penalty u(x), entering such region yields a relatively high cost. L(x, v) = 1 2v 2+u(x), u(x) = 100. 6 37
and nonlocal interactions, uneven splitting caused by asymmetric nonlocal kernel (nonpotential). I Density constraint over space ⇥ time: ⇢(x, t) cUpperBound, (x, t) 2 ⌦ ⇥ [0, T]. I Obstacles (yellow) modeled via density constraint:⇢(x, t) = 0, x 2 ⌦ Obstacles. 7 38
minimize ⇢,v Z T 0 Z Td 1 2⇢vT GM vdx + F(⇢(·, t)) dt + G(⇢(·, T)) subject to ⇢t + r · (⇢v) = 0 ⇢(·, 0) = ⇢0, (1) We call (1) the forward problem and call the problem of recovering the ground metric GM and the cost functionals F from the observations of ⇢ and v the inverse problem. Particularly, in maximum mean discrepancy (MMD) case, F(⇢(·, t)) = Z Td 1 2⇢(x, t)(K ⇤ ⇢)(x, t)dx. (2) We are interested in the convolution kernel K in (2). 2 40
the game I v: velocity ﬁeld of the agents Target parameters (parameters to be recreated): I Ground metric GM : a metric function on the sample space, depicting the geometric contour of the transport space. I convolution kernel in MMD K: depicts the pairwise impact between the particles, which relies on nothing but the relative distance of any two players. 3 41
↵ 2 k⇢ ˆ ⇢k2 + 2 kv ˆ vk2 + ↵0 2 k⇢0 ˆ ⇢0 k2 + k⇢T ˆ ⇢T k2 + p kr✓kp p subject to r · (⇢v) + ⇢t = 0 (GM v)t + r ✓ 1 2 vT GM v ⇢ F(⇢) ◆ = 0 @(GM v)i @xj = @(GM v)j @xi , i 6= j Z si(ˆ xi,·) (GM v)idSx = 0, ˆ x i 2 Td 1 , i = 1, 2, . . . , d. (3) The objective function consists of L2 distance between density distribution, velocity ﬁeld and their observations, and an additional Hp norm on target parameter ✓. The constraints restrict the solution to satisfy MFG PDEs. 4 42
with primal-dual algorithm and Bregman iteration. Numeric result by applying primal-dual algorithm: 0 50 0.5 50 1 0 0 0 50 0.5 50 1 0 0 0 50 0.1 50 0.2 0 0 Figure: Ground metric recreated in 2-dimensional basis under H 2 regularization. Parameters in the model were taken as ↵ = 1000/kˆ ⇢k2 , = 1/kˆ vk2 , ↵0 = 0, = 10 2. From left to right: g0 learned from the observations, real ¯ g0 , the absolute di↵erence between the learned g0 and the real data |g0 ¯ g0 |. 5 43
10 0 0 0 20 0.5 20 10 1 10 0 0 0 20 0.1 20 10 0.2 10 0 0 Figure: Convolution kernel recreated in 2-dimensional basis under H 2 regularization. Parameters in the model were taken as ↵ = 100/kˆ ⇢k2 , = 1/kˆ vk2 , ↵0 = 0, = 10 3 From left to right: kernel K(x, 0) learned from the observations, real kernel Kr(x, 0), the absolute di↵erence between the learned parameter and the real one |K(x, 0) Kr(x, 0)|. 6 44
model (Kermack and McKendrick, 1927) 8 > > > > > < > > > > > : dS dt = SI dI dt = SI I dR dt = I where S, I,R : [0, T] ! [0, 1] represent the density of the susceptible population, infected population, and recovered population, respectively, given time t. The nonnegative constants and represent the rates of susceptible becoming infected and infected becoming recovered. 5 49
,the spatial SIR model is considered: 8 > > > > > > > < > > > > > > > : @ t ⇢ S (t, x) + ⇢ S (t, x) Z ⌦ K(x, y)⇢ I (t, y)dy ⌘2 S 2 ⇢ S (t, x) = 0 @ t ⇢ I (t, x) ⇢ I (x) Z ⌦ K(x, y)⇢ S (t, y)dy + ⇢ I (t, x) ⌘2 I 2 ⇢ I (t, x) = 0 @ t ⇢ R (t, x) ⇢ I (t, x) ⌘2 R 2 ⇢ R (t, x) = 0 Here ⌦ is a given spatial domain and K(x, y) is a symmetric positive deﬁnite kernel modeling the physical distancing. E.g. R Kd⇢ I is the exposure to infectious agents. 6 50
. Deﬁne the Lagrangian functional for Mean ﬁeld game SIR problem by L((⇢ i , m i , i ) i=S,I,R ) =P(⇢ i , m i ) i=S,I,R Z T 0 Z ⌦ X i=S,I,R i ✓ @ t ⇢ i + r · m i ⌘2 i 2 ⇢ i ◆ dxdt + Z T 0 Z ⌦ I ⇢ I K ⇤ ⇢ S S ⇢ S K ⇤ ⇢ I + ⇢ I ( R I )dxdt. Using this Lagrangian functional, we convert the minimization problem into a saddle problem. inf (⇢i,mi)i=S,I,R sup ( i)i=S,I,R L((⇢ i , m i , i ) i=S,I,R ). 16 52
⇢ i (0, ·) (i = S, I, R) Output: ⇢ i , m i , i (i = S, I, R) for x 2 ⌦, t 2 [0, T] While relative error > tolerance ⇢(k+1) i = argmin ⇢ L(⇢, m(k) i , (k) i ) + 1 2⌧i k⇢ ⇢(k) i k2 L2 m(k+1) i = argmin m L(⇢(k+1), m, (k) i ) + 1 2⌧i km m(k) i k2 L2 (k+ 1 2 ) i = argmax L(⇢(k+1), m(k+1) i , ) 1 2 i k (k) i k2 H2 (k+1) i = 2 (k+ 1 2 ) i (k) i end 17 53
status of populations, going beyond S, I, R. I Construct discrete spatial domain model, including airport, train station, hospital, school etc. I Propose inverse mean ﬁeld SIR problems. Learning parameters in the model by daily life data. I Combine mean ﬁeld game SIR models with AI and machine learning algorithms, including APAC, Neural variational ODE, Neural Fokker-Planck equations, etc. 22 57