Self-learning Monte Carlo method andall optical neural network
Junwei Liu (劉軍偉)
Department of Physics
Hong Kong University of Science and Technology
Deep learning and Physics, YITP, Kyoto, Japan, 2019

The general motivation•

As well known, there are great developments in machinelearning techniques such as principle component analysis,
deep neural networks, convolutional neural networks,
generative neural networks, reinforcement learning and so on.

•By using these methods, there are also many great
achievements such as image recognition, auto-pilot cars andAlpha-GO.

•Can we use these methods in physics to solve some problems?
If so, what kind of problems can we solve and how?

2/46

Part 1: Self-learning Monte Carlo (ML Physics)

Monte Carlosimulation
3/46

Collaborators and References

Liang Fu

Yang Qi(MIT)(Fudan)Ziyang Meng Xiaoyan Xu(IOP)(UCSD)

Yuki Nagai(JAEA)

Huitao Shen(MIT)

References
1.
2.
3.
4.

5.

Self-Learning Monte Carlo Method, PRB 95, 041101(R) (2017)
Self-Learning Monte Carlo Method in Fermion Systems, PRB 95, 241104(R) (2017)
Self-Learning Quantum Monte Carlo Method in Interacting Fermion Systems, PRB
96, 041119(R) (2017)
Self-learning Monte Carlo method: Continuous-time algorithm, PRB 96, 161102(2017)
Self-learning Monte Carlo with Deep Neural Networks, PRB 97, 205140 (2018)

Related work from Prof Lei Wang’s group at IOP
1.
2.
4/46

Accelerated Monte Carlo simulations with restricted Boltzmann machines, PRB 95,035105 (2017)
Recommender engine for continuous-time quantum Monte Carlo methods, PRE 95,031301(R) (2017)

Monte Carlo methods•Consider a statistical mechanics problem
𝑍 = ෍ 𝑒 −𝛽𝐻[𝐶] = ෍ 𝑊(𝐶)𝐶𝐶𝑂 = ෍ 𝑂 𝐶 𝑒 −𝛽𝐻 𝐶 ൗ𝑍 = ෍ 𝑂 𝐶 𝑊(𝐶)ൗ𝑍

•

𝐶𝐶

Pick up 𝑁 configurations (samples) in the configuration space
𝐶 based on the importance 𝑊 𝐶𝑖 /𝑍. Then we can estimate
observables 𝑂 as𝑁𝑂 ≈ ෍ 𝑂 𝐶𝑖 ൘𝑁𝑖=1

•

5/46

The statistical error is proportional to 1Τ 𝑁. In high dimension(d>8), Monte Carlo is the most important and sometimes the
only available method to perform the integral/summation.

Quantum Monte Carlo methods•Consider a quantum system characterized by 𝐻𝑍 = ෍ 𝜓 𝑒 −𝛽𝐻 𝜓𝜓Method 1: Trotter decomposition
𝑍 = ෍ 𝜓 𝑒 −𝛽𝐻 𝜓𝜓=σ𝜓𝜓1 𝜓2…𝜓𝑀 𝜓1 𝑒 −𝛽𝐻/𝑀 𝜓2 𝜓2 𝑒 −𝛽𝐻/𝑀 𝜓3 … 𝜓𝑀 𝑒 −𝛽𝐻/𝑀 𝜓1
Method 2: serial expansion
𝑍 = σ𝜓 𝜓 𝑒 −𝛽𝐻 𝜓 =σ𝜓 𝜓 σ𝑛 −𝛽𝐻 𝑛 /𝑛! 𝜓• Map the N-dimensional quantum model to be a (N+1)dimensional “classical” model, and then use Monte Carlo
method to simulate this “classical” model.
6/46

Markov chain Monte Carlo (MCMC)•MCMC is a way to do important sampling based on thedistribution 𝑊(𝐶). Configurations are generated one by one bythe following approach:
1. first propose the next trial configuration 𝐶𝑡
2. If 𝑅𝑎𝑛𝑑() ≤ 𝑻(𝑪𝒊 → 𝑪𝒕 ), then the next configuration 𝐶𝑖+1 = 𝐶𝑡
Otherwise, 𝐶𝑖+1 = 𝐶𝑖
3. Repeating step 1 and step 2⋯ → 𝐶𝑖−1 → 𝐶𝑖 → 𝐶𝑖+1 → ⋯•7/46

It is clearly that the next step only depends on the last stepand transition matrix 𝑇, and we can demonstrate that thedetailed balance principle guarantees that the Markov
process will converge to the desired distribution
𝑇(𝐶 → 𝐷) 𝑊 𝐷=𝑇(𝐷 → 𝐶) 𝑊 𝐶
N. Metropolis, J. Chem. Phys. 21, 1087 (1953)

Metropolis-Hastings algorithm•

In Metropolis Algorithm, transition matrix is chosen as
𝑊 𝐷𝑇 𝐶 → 𝐷 = min 1,𝑊 𝐶
N. Metropolis, et. al J. Chem. Phys. 21, 1087 (1953)•The transition probability can be further written as
𝑇 𝐶→𝐷 =𝑆 𝐶→𝐷 𝑝 𝐶→𝐷
𝑆 𝐶 → 𝐷 :proposal probability of conf. 𝐶 from conf. 𝐷
𝑝 𝐶 → 𝐷 :acceptance probability of configuration 𝐶•

In Metropolis-Hastings Algorithm, the acceptance ratio is𝑊 𝐷 𝑆 𝐷→𝐶𝑝 𝐶 → 𝐷 = min 1,𝑊 𝐶 𝑆 𝐶→𝐷8/46

W. H. Hastings, Biometrika 57, 97 (1970)

Independent Samples•

As well known, for thestatistic measurements,
only the independent
samples matter.

•However, since the configurations are generated sequentially
one by one as random walk in configuration space, it isinevitable that the configurations in the Markov chain arecorrelated with each other and not independent.

•The configurations generated by different Monte Carlo
methods have different autocorrelations.

9/46

Different update algorithms•

Too small step length: small difference, high acceptance.

•Too large step length: big difference, low acceptance.
𝒆−𝜷(𝑬 𝑫 −𝑬(𝑪))•10/46

Ideal step length: big difference and high acceptance,
exploring the low-energy configurations.

How to justify different Monte Carlo methods?
Time consumption 𝑡𝑙𝑐 to get two statistically independentconfigurations• Auto-correlation time 𝑙𝑐 : the number of steps to get
independent configurations
 Bigger differences  independent, but low acceptance ratio
 Similar energies  high acceptance ratio, but notindependent• Time consumption 𝑡 to get one configuration (mainly for thecalculation of weight)

Self-learning Monte Carlo methods are designed to improve
both parts, thus can speed up the calculations dramatically.
11/46

Local update• Local Update

𝑆 𝐶→𝐷 =𝑆 𝐷→𝐶

1=𝑁• Acceptance ratio

𝑊 𝐷 𝑆 𝐷→𝐶
𝛼 𝐶 → 𝐷 = min 1,𝑊 𝐶 𝑆 𝐶→𝐷= min 1, 𝑒 −𝛽(𝐸 𝐷 −𝐸(𝐶))• Very general: applies to any model
12/46

N. Metropolis, et al J. Chem. Phys. 21, 1087 (1953)

Critical slowing down•

Dynamical relaxation time diverges at the critical point:convergence is very slow in the critical system.

•For 2D Ising model, autocorrelation time 𝜏 ∝ 𝐿𝑧 , 𝑧 = 2.125

13/46

How to get high acceptance ratio?

𝑊 𝐷 𝑆 𝐷→𝐶
𝛼 𝐶 → 𝐷 = min 1,𝑊 𝐶 𝑆 𝐶→𝐷𝑊 𝐷
𝑆 𝐶→𝐷=𝑊 𝐶
𝑆 𝐷→𝐶

14/46

Global update -- Wolff algorithm in Ising model

1. Randomly choose one site 𝑖
2. If the adjacent sites have the same status, then add them inthe cluster 𝐶 with the probability 𝑆 𝑖 → 𝑗 = 1 − 𝑒 −2𝛽𝐽
3. Repeat step 2 for all the sites in the cluster 𝐶
4. Change the status of all the sites in the cluster 𝐶
15/46

Swendsen and Wang, Phys. Rev. Lett. 58, 86 (1987)
U. Wolff, Phys. Rev. Lett. 62, 361 (1989)

Reduce critical slowing down

16/46

Swendsen and Wang, Phys. Rev. Lett. 58, 86 (1987)

Local update and global update

Local update

Global update• Locally update configuration • Globally update configurationby changing one site per MC
by simultaneously changingstep
many sites per MC step
• Very general• High efficiency• Inefficient around the phasetransition points (criticalslowing down)• Designed for specificmodels and hard to begeneralized to other models

17/46

How to get high acceptance ratio?

𝑊 𝐷 𝑆 𝐷→𝐶
𝛼 𝐶 → 𝐷 = min 1,
𝑊 𝐶 𝑆 𝐶→𝐷𝑊 𝐷
𝑆 𝐶→𝐷=𝑊 𝐶
𝑆 𝐷→𝐶𝑊 𝐷
𝑆 𝐶→𝐷
≈𝑊 𝐶
𝑆 𝐷→𝐶
18/46

VeryHard

My initial naïve idea••Use machine learning to learn some common features forthese “important” configurations, and then generate newconfigurations based on these learned features.
It seems to be good, but it does not work, because we don’tknow•
19/46

𝑆 𝐶→𝐷, and we cannot calculate the right acceptance
𝑆 𝐷→𝐶probability.
However, this idea tells us that there are other important
information in the generated configurations, besides 𝑂 ≈σ𝑁𝑖=1 𝑂 𝐶𝑖 Τ𝑁 based on them.

Hidden information in generated configurations•

The generated configurations have a similar distribution close
to the original distribution !!!

•It is too obvious and seems to be too trivial, but we do not use
this seriously besides calculating the average of operators.

•

Can we further use it and how?

•The answer is YES and we can do it in self-learning Monte
Carlo method.

20/46The right way to use the hidden information in
generated configurations

Yang Qi(Fudan)

Critical properties

Learn

Guide
Universality class

Effective Model
21/46

Core ideas of self-learning Monte Carlo• Learn an approximated simpler model Having efficient global update methods
 Evaluating the weights faster

First learn

How totrain?

• Use the simpler model to guide us tosimulate the original hard model

Then earn

How topropose?
22/46

J. Liu, et al PRB 95,041101(R) (2017)

SLMC in Boson system
•

23/46

Original Hamiltonian has both two-body and four-bodyinteractions, and we do not have global update method.

Fit the parameters in 𝑯eff•

Effective Hamiltonian••Generate configurations with local update at T=5>Tc, awayfrom the critical points
Perform linear regression to fit the parameters in 𝐻eff•

Generate configurations with reinforced learning at Tc

24/46

Build the cluster based on 𝑯eff•

Self-learning update: cluster is constructed using Wolff update
obeying detailed balance principle of effective model 𝑯eff
𝑊eff 𝐶𝑆 𝐷→𝐶=𝑆 𝐶→𝐷𝑊eff (𝐷)• Acceptance ratio:
𝑊 𝐷 𝑊eff 𝐶𝑊 𝐷 𝑆 𝐷→𝐶
𝛼 𝐶 → 𝐷 = min 1,= min 1,𝑊 𝐶 𝑆 𝐶→𝐷
𝑊 𝐶 𝑊eff (𝐷)= min 1, 𝑒

25/46

−𝛽 𝐸 𝐷 −𝐸

eff 𝐷 − 𝐸 𝐶 −𝐸eff 𝐶

For different sizes 𝑳 × 𝑳••
26/46

SLMC is 20~40 fold faster than the local update.
We can easily get the result for 320 × 320.
J. Liu, et al PRB 95, 041101(R) (2017)

How to justify different Monte Carlo methods?
Time consumption 𝑡𝑙𝑐 to get two statistically independentconfigurations• Auto-correlation time 𝑙𝑐 : the number of steps to get
independent configurations
 Bigger differences  independent, but low acceptance ratio
 Similar energies  high acceptance ratio, but notindependent• Time consumption 𝑡 to get one configuration (mainly for thecalculation of weight)

Self-learning Monte Carlo methods are designed to improve
either or both parts depending on models, thus can speed
up the calculations dramatically.
27/46

Weight (Free energy) in Fermion systems• In general fermion systems, the partition function can becalculated based on the following form

• We have to deal with heavy matrix operations to calculate theweight (free energy), which is very time-consuming 𝑂(𝑁 3 ).

• However, the weight only depends on the configuration ofbosonic field, which can be described by pure Boson models.
28/46

Complexity of SLMC

• Complexity to obtain two uncorrelated configurations 𝑁𝜏 × 𝑁• Complexity to obtain two uncorrelated configurations
29/46

3 = 𝜏𝑁 4𝑁𝜏 + 𝑁 3 /𝑝

Example: double exchange model• Magnetic order in metal• Free Fermions are coupled to classical Heisenberg spins• Weight• The 𝐸𝑛 are obtained by ED of 𝐻 for given configurations 𝜙• Computational complexity in conventional method: 𝑂(𝜏𝐿𝑑 𝐿3𝑑 )30/46

Learning the effective model• Effective model• Reproduce the RKKY feature• Reproduce the right Boltzmann distribution
31/46

J. Liu, et al PRB 95, 241104(R) (2017)

No critical slowing down

• Autocorrelation time is much shorter in SLMC than conventional
method• There is no critical slowing down in SLMC
• Computational complexity in conventional method: 𝑂(𝜏𝐿𝑑 𝐿3𝑑 )• Computational complexity in SLMC: 𝑂 𝑙𝑐 ≈ 𝜏𝐿𝑑 + 𝑂(𝐿3𝑑 )• SLMC reduces the complexity by 𝑂(𝜏𝐿𝑑 ), and can be more than
1000 times faster
32/46

J. Liu, et al PRB 95, 241104(R) (2017)

Example II: SLMC in DQMC (Interacting Fermion)• In conventional DQMC, the computational
complexity in conventional method:𝑂(𝜏𝛽𝑁 3 )Ziyang Meng(IOP)• In SLMC with cumulative update, thecomplexity is 𝑂(𝛾𝛽𝑁𝑙𝑐 + 𝑁 3 + 𝛽𝑁 2 )a. Cumulative update 𝑂(𝛾𝛽𝑁𝑙𝑐 )b. Detail balance 𝑂(𝑁 3 )c. Sweep Green’s function 𝑂(𝛽𝑁 2 )33/46

Xiaoyan Xu(UCSD)

X. Y. Xu, et al PRB 96, 041119(R) (2017)

Results for a particular model

For the first time, we canachieve large system size as
100×100 in two dimension

34/46

X. Y. Xu, et al PRB 96, 041119(R) (2017)

Example III: SLMC in CT-AUX
Single impurity Anderson model

Yuki Nagai(JAEA)For the single impurity Anderson model
in low temperature or large interaction,
one can achieve speedup𝑛𝑂≈ 𝑂 𝛽𝑈𝑚𝑐 𝜏𝑆𝐿

35/46

Yuki Nagai, et al PRB 96, 161102 (2017)

SLMC with Deep Neural Networks

NeuralNetworks

NeuralNetworks

36/46

H. Shen, J. Liu, L. Fu, PRB 97, 205140 (2018)

Asymmetric Anderson model with a single
impurity• 𝑑መ 𝜎 and 𝑐𝑘Ƹ are the fermion
annihilation operator for theimpurity and for theconduction electrons, and𝑛ො 𝑑,𝜎 ≡ 𝑑መ 𝜎+ 𝑑መ 𝜎 , 𝑛ො 𝑑 = 𝑛ො 𝑑↑ + 𝑛ො 𝑑↓• By Hubbard-Stratonovich transformation, we can decouple theinteraction parts as an auxiliary bosonic field coupled to freefermions• Then we can use Monte Carlo to sample the auxiliary field
configuration space, thus simulating the original asymmetric
Anderson model.

• It is not easy to write an explicit effective model to describe thebosonic action here, but we use neural networks to do it.
37/46

Fully connected neural network

38/46

Insights from fully connected NN• Translation symmetry• Much less parameters• Easy to extend

39/46

H. Shen, J. Liu, L. Fu, PRB 97, 205140 (2018)

Results of convolutional network

Huitao Shen(MIT)• DNNs can generally
work for different
chemical potential
and temperature.

H. Shen, J. Liu, L. Fu,
PRB 97, 205140 (2018)
40/46

Part 2: All optical neural network (Physics ML)
Shengwang Du(HKUST)

41/46• At the speed of light• Intrinsic infinite parallel calculations• Robust against local errors

Spatial light modulator: linear transformation• Use grating to control therelation between the input
light and the output light linear transformation
𝐼𝑜𝑢𝑡 = 𝑊𝐼𝑖𝑛
42/46

Test the linear transformation
L2

L3
FM

SLM1

SLM2

L1

Camera1

Camera2

Coupling

43/46

L4

M

Non-linear function in optics

Electromagneticallyinduced transparency

44/46

Test the fully connected neural network• We can reproduce well the results from the fully connectedneural networks (16, 4, 2) in computer
45/46

Ying Zuo, et al. Optica 6, 1132 (2019)

Summary
1.

SLMC can be generally applied in both bosonic and fermionicsystems to speed up MC simulations.

2.

SLMC is a bridge between the numerical simulation and analytic
studies for many-body physics:•

Better understanding from analytic studies can give us better
effective model and better efficiency.

•The learned effective Model in SLMC is a good starting point,and can give us more insights for analytic studies.

3.

SLMC also offers us a framework to integrate machine learningtechniques into Monte Carlo simulations.

4.

We have already realized the AONN in experiments.

46/46