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
@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
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 ﬂow is a gradient ﬂow of entropy in optimal transport metric. 3 8
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
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
inner product g(u): TuM ⇥ TuM ! R is given below. For any 1 , 2 2 TuM, deﬁne 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
: 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 inﬁmum is taken among all density functions u: [0, 1] ⇥ ⌦ ! R, vector ﬁelds 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 ﬁxed initial and terminal density functions u0 , u1 2 M(⌦). 8 13
ﬂows) 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-ﬁeld information metric forms g(u)( 1, 2) = Z ⌦ (r 1(x), r 2(x))u(x)dx. Here the Hamiltonian ﬂow in (M(⌦), g) satisﬁes 8 < : @tu + r · (ur ) = 0 @t + 1 2kr k 2 + uF(u) = 0. 10 15
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 speciﬁc norm H which allows step sizes to be independent of an operator A. ⌧kAk2 < 1 (classical PDHG) ! ⌧ < 1 (G-Prox PDHG). 12 17
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
scalar conservation law deﬁned 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
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 ﬁnd 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
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 deﬁned 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
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 ﬁrst-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
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 0lNt 1 1jNx 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
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
= @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
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
equation of the following form: 8 > < > : @t⇢ + @xm = 0, @tm + @x( m2 ⇢ ) + @xP(⇢) = @x(µ(⇢)@x m ⇢ ), Here, ⇢: density; m := ⇢v: ﬂx; 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
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
Equilibrium on ﬂat 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-ﬁeld games to manifolds, and show that with operators on manifolds, the Nash Equilibrium satisﬁes 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
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
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
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
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
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