1

ベイズ推定を用いた画像パター
ンからの支配方程式の選択
義永

那津人

WPI-AIMR, 東北大学
MathAM-OIL, 産総研

共同研究者：
徳田 悟（九州大学）

2

Machine Learning, Data-Driven Science, Inverse
Problem, Optimisation• forward problem
input

model

structuretemperaturefunction,
equationx

outputyf ( x)elasticity
structure (pattern)• inverse problem
datax, ymodelestimationy  f ( x)3

We do not know the true model…
• Parametrised modelf ( x)   0  1 x   2 x 2  3 x3 ▫ polynomialsf ( x)  1 cos x   2 cos 2 x  ▫ Fourier expansion▫ deep learning
f w,b ( x)   Wn ... (W1 x  b1 )   bn 

• Learning▫ supervised▫ unsupervised{x, y}  y  f ( x)x  f ( x) such that x is linearly

discriminable・normal dist. to Boltzmann dist.
(z)
p(x)
▫ generative
・deep learninginterpretability or predictability

4Contents• Discovery of a governing equation from timeseries data (Nathan Kutz)• Inverse structural design for molecular
dynamics and Monte Carlo simulations(Thomas Truskett)• Bayesian modelling of partial differential
equations from one snapshot of pattern(arxiv:2006.06125)

5

Discovery of a Model (Brunton-Proctor-Kutz (2016))• Ordinary differential equationsx (t )  f ( x(t ))x (t )  f ( x(t ))x(t )   0  1 x   2 x 2  3 x 3  • Partial differential equations

 t u ( x, t )  f (u ( x, t ))

 t u  f (u ,  x u ,) t u   0  1,0u   2,0u  1,1 x u  1,2  xx u 2Can we estimate a model from data?
Parametrise r.h.s. by  and perform regression

6Sparse Identification of Dynamical systems (SINDy,PDE-FIND)
model

 t u   0  1,0u   2,0u  1,1 x u  1,2  xx u 2data

{xi , ti , ui ( xi , ti )}i 1,N
time derivatives

sparse regression

 t U  (U)  candidate terms2  t u1   u1
 
u2t   u2
    u3
 u t N   x u1
 x u2
 x u3u1    1,0  
u22    1,1 
u32     2,0      parameters  arg min    t U    (U)    2sparsity

7

Example: reaction-diffusion equations
ground-truth model

Rudy-Brunton-Proctor-Kutz (2017) t u ( x, t )  u  (u 2  v 2 )u  (u 2  v 2 )v  0.1uu

 t v( x, t )  v  (u 2  v 2 )v  (u 2  v 2 )u  0.1vv

sparsityestimation
 t u ( x, t )  1.00u  1.00u 3  1.00v3  1.00uv 2  1.00u 2 v  0.100 xx u  0.100 yy u
 t v( x, t )  1.00v  1.00v 3  1.00u 3  1.00uv 2  1.00u 2 v  0.100 xx v  0.100 yy v

 u2  
 uv 
0  
 x u   

8Discovering a Model

Lorentz modelx    y  x 

Korteweg–de Vries (KdV) equationy  x    z   y t u  6u x u   xxxu  0z  xy   zBurgers equationx

 t u  u x u   xx u

space

Nonlinear Schrödinger
12i t   xx    0timeu (x, t ) v(x, t ) ( x, t )

time

2

Kuramoto–Sivashinsky equation12
 t u    xu    xx u   xxxx u  0
2reaction-diffusion equation
Navier-Stokes equation
Re   t   v      zxyBrunton-Proctor-Kutz (2016)
Rudy-Brunton-Proctor-Kutz(2017)

9Remarks• Equation discovery for time series has a long history▫ system identification (Crutchfield-McNamara 1987, Bär et al. 1999,Muller-Timmer 2004)▫ Finding an equation of motion for a double-pendulum from data
(Bongard-Lipson 2007, Schmidt-Lipson 2009)▫ Sparse regression is simple and reasonably powerful, and looks
widely applicable.

• Application▫ time-dependent parameters (group sparsity)▫ hidden variables• Not always work

▫ Precise evaluation of time and spatial derivatives is necessary.
▫ Not robust against noise.
▫ Need a lot of data (enough information including relaxation process).
▫ What you get is a normal form, not necessarily a ground-truthequation.
▫ Assume the state noise, not observation noise10

Inverse Structural Design• forward problem▫ molecular dynamics simulations,dynamics, Monte Carlomodel

initial x (0)condition i

Brownianstationarystate p ({xi })xi  U ({xi })  i
Brownian dynamics+ potential energy• inverse problem
target structurexi ~ ptgt ({xi })modelestimationU ({xi })ptgt ({xi })  p({xi }|  )11Inverse Structural Design
partition function(normalisation)

Boltzmann distributionp ({xi }|  ) e  U ({ xi })    e  U ({ xi }) dxi

i

Kullback-Leibler (KL) divergence
DKL  ptgt ({xi }) || p ({xi }|  )    ptgt ({xi }) lnptgt ({xi })p({xi }|  )dxi

Find  satisfying  DKL  ptgt ({xi }) || p({xi }|  )   0
 ( k 1)   ( k )    ln p ({xi }|  ) P ({ x }) tgti  ( k ) U ({xi }|  ) P ({ x }| )   U ({xi }|  ) P ({ x })
i

tgti

12

Inverse Structural Design
Compare radial distribution functions
( k 1)

targetstructureg tgt (r )(k )   dr  g (r |  )  g tgt (r )   u (r |  )  ( k )
0
 
Monte
generated
Carlo
model
structure
U ({x })g (r |  )iestimation

Lindquist-JadrichTruskett (2018)

13Remarks• Parameter estimation▫ Nothing but parametrised potential estimation▫ Boltzmann machine with a specific form of thepotential function.

• Sampling▫ Generating samples from the probability distribution
is not easy.
▫ Sophisticated Monte Carlo methods (annealing)• Quality of the estimation▫ Point estimate; most likely potential is estimated.
▫ No uncertainty quantification▫ Stable to the change of parameters?

14

Bayesian modelling of partial
differential equations (BM-PDE)• forward problem▫ Stationary patterndifferential equations
initial  0 ( x)condition

generated

PDE model t (x, t )  f  (x, t ) bypartial

stationarypattern  s ( x)• inverse problem
target pattern * ( x)PDE modelestimation

 t (x, t )  f   (x, t )  * (x)   s ,  (x) up to symmetry transformation

15

Pattern Formation
block copolymers

double gyroid

phenomenological model
partial differential equations

wide application




t

free energy functional


2
1

      f GL    (x)  G (x, x ')  (x ')  
 22 '











· reaction-diffusion· fluid· phase-field model· crystals(phase-field crystal)16

Order Parameters• Requirement▫ Patterns are invariant undertranslation and rotation.
▫ Patterns respect their symmetry.
▫ Small noise must be eliminated.

• Order Parameter▫ Fourier transformation  k    (x) ▫ Expand the coefficientsAl    k eilk2Dk
2D invariant  l  Al

 cos  k k k

sin
k 

17

Can we estimate the best model from asnapshot?
• Very limited information▫ No information about time series▫ Symmetry is not known.
▫ Not unique; Regularisation is necessary.
▫ Quantify quality of the estimation; we don’t know theground-truth.
2D quasi crystalpattern ( x)*?

model   (x, t ) (x, t )t parameters

NY-Tokuda (2020)

18

Bayesian modelling of partial
differential equations (BM-PDE)

・Make a family of PDE models (the number of length scales)
・Bayesian modelling + sampling by replica exchange Monte Carlo method
Parameter estimation from the posterior distribution, and model selection by
marginal likelihood

19

Phase field crystal model
Nonlinear partial differential equations(gradient system with mass conservation)

 t  

For quasicrystals
Lifshitz-Petrich (1997)Swift-Hohenberg-type free energy
21

1 2
 

[ ]      2      1   q 2    3   4 

234 
 2





second lengthscale
dodecagonal patternq  2 cos( /12)

quasicrystal

N  128

Elder, et al (2002)

20

Family of models• Family of linear operator▫ Parametrised by characteristic wave numbers andtheir eigenvalues ( m )       t k

 ,kkk2
a0 2 2s0 2 2( m1 )22k  2q0  4 k q0  k
 , k k    4 k
q0q0



one length scale

two length scalethree length scale

21

Bayes theorem
likelihood  prior distribution
posterior distribution 
normalisation constanttarget patternparameter estimation

p   |  ,  , m *

generated patternparametersp  * |   ,  , m p (  | m)p * |  , mnoise level
model selection model(inverse temperature)observation noise *    noise

uncertainty quantification(best noise level anduncertainty of parameters)

22Generating a pattern through the estimation
target pattern(real space)target pattern(Fourier space)

generated pattern(estimated modeland parameters)

12 spots
dodecagonal quasicrystal

23

Model Selection and Parameter Estimation
model selection from log marginal
likelihood→ two-length scales are necessary
parameter estimation fromposterior distribution
ground truth

24

Three-Dimensional Patterns
double gyroid generated from a function

pattern generated from estimated model and parameters

25

Various patterns in 2D and 3D

NY-S.Tokuda (2020)

26Summary• Discovery of a governing equation▫ Not difficult than it looks▫ Becoming feasible▫ Yet, we need quantification of its quality• Bayesian modelling of partial differential equations▫ BM-PDE works well for various patterns.
▫ It works not only for a target pattern with the groundtruth, but also for a pattern without it.
▫ It is also robust against noise (~20%).
Natsuhiko Yoshinaga and Satoru Tokuda“Bayesian modelling of partial differential equations from one
snapshot of pattern“, arxiv:2006.06125