2023-07-24 Extreme Universe Collaboration - 18th Colloquium

Markov-Chain Monte Carloin Tensor-Network Representation
Collaboration

Extreme Universe
The 18th COLLOQUIUM

July 24 th (Mon.) ONLINETALK 10:00 - 11:00 (JST)July 24 th (Mon.) 1:00 - 2:00 (UTC)
July 23 th (Sun.) 21:00 - 22:00 (EDT)ONLINE COFFEE TIME
11:00 - 12:00 (JST)Registration required (click HERE)
Extreme Universe, JAPANSpeakerProf. Synge Todo
The University of Tokyo

Title

Markov-Chain Monte Carloin Tensor-Network Representation

Synge Todo / 藤堂眞治 <wistaria@phys.s.u-tokyo.ac.jp>
Department of Physics, University of Tokyo

2023

Abstract
The partition functions of various classical and quantum lattice models can be represented as
tensor networks. However, the exact contraction of a tensor network is generally exponentiallyexpensive, and some approximation, such as low-rank approximation based on singularvalue decomposition, is usually required. Recently, we proposed a new tensor contraction
method based on Monte Carlo sampling. The proposed method combines the stochastic
basis transformation of tensors with the Markov chain Monte Carlo framework.
It can entirely remove the systematic error due to a finite bond dimension of the low-rankapproximation while keeping the high accuracy of the tensor-network method. We also
demonstrate how the proposed method works for systems with negative (or complex)
weights, where the standard Markov chain Monte Carlo suffers from a severe sign problem.

MEXT -KAKENHI- Grant-in-Aid for Transformative Research Areas (A)The Natural Laws of Extreme Universe -A New Paradigm for Spacetime and Matter from Quantum Information-

2

Agenda• Introduction• Markov-chain Monte Carlo and negative sign problem• Tensor network representation and approximate contraction• Sampling approach in tensor network contraction• Projector formulation of tensor network method• Sampling projectors• Monte Carlo in tensor network representation• Sequential Monte Carlo in tensor network representation• Markov-chain Monte Carlo in tensor network representation• Results• Summary3

Stochastic approach in computational physics• “Our brains are just not wired to do probability problems very well”- Persi Diaconis• Markov chain Monte Carlo - Metropolis et al (1953)• can sample from an arbitrary probability distribution• perfect sampling (“coupling from the past") - Propp and Wilson (1996)• samples perfectly independent samples from Markov chain• extended ensemble method - Hukushima and Nemoto (1996), Wang andLandau (2001)• realizes equilibrium immediately after quench• samples extremely rare events ( ∼ 10−100)• O(N ) method for long-range interacting system - Fukui and Todo (2009)• exact O(N ) sampling (and energy measurement) instead of O(N 2)• MCMC without detailed balance - Suwa and Todo (2010), Michel et al(2014)• diffusive dynamics → ballistic dynamics

4

Advances in Markov chain Monte Carlo• Representation (de nition of “con gurations” and “weighs”)• path integral representation for quantum Monte Carlo (1976), Bayesianinference (1990)...

• Choice of ensemble• extended ensemble method: multicanonical MC (1991, 2001),
exchange MC (1996), lifting (2000)...

• Generation of set of candidate con gurations• non-local (cluster) updates: Swendsen-Wang (1987),
Hamiltonian MC (1987), loop (1993), worm (1998)...

• Choice of transition kernel (probabilities)• Metropolis, heat bath (Gibbs sampler), over-relaxation (1987),
irreversible kernel (2010), event-chain (2013)...

• Algorithm for generating a con guration according to transition
probabilities

fi

fi

fi

fi

• N-fold way (rejection free) (1975), Walker’s method (1977, 2019),
order-N algorithm (1995, 2009)...

5

Negative sign problem• In the path-integral representation for frustrated magnets, fermionicsystems, real-time dynamics, the sample weights become negative(or even complex)• average sign becomes exponentially small by cancellation at lower
temperatures, longer time, and/or larger system sizes• speci c heat of antiferromagnetic Heisenberg model on kagomé lattice
1.4

speci c heat

1.2
1

average sign

0.8
0.6
0.4
0.2

fi

fi

0

T
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

6

Unitary dynamics in quantum circuits

| 0⟩

H

| 0⟩
| 0⟩

H

H

1
1
|
|
|
|
|
|
000⟩
⇒
(
0⟩
+
1⟩)
0⟩(
0⟩
+
1⟩)
=
( | 000⟩ + | 001⟩ + | 100⟩ + | 101⟩)
•
2
2
1
1
⇒ ( | 000⟩ + | 001⟩ + | 110⟩ + | 111⟩) ⇒ ( | 000⟩ + | 001⟩ + | 110⟩ − | 111⟩)
2
2
1⇒( | 00⟩ + | 11⟩)( | 0⟩ + | 1⟩)+( | 00⟩ − | 11⟩)( | 0⟩ − | 1⟩)
2 2
1=( | 000⟩ + | 001⟩ + | 110⟩ + | 111⟩ + | 000⟩ − | 001⟩ − | 110⟩ + | 111⟩)
2 2
1=( | 000⟩ + | 111⟩)
2

7

Path sampling of quantum circuits

| 0⟩

H

| 0⟩
| 0⟩

H

| 000⟩• State transition diagram
000

H

000

100

000

000

000

000

001

001

001

001

110

110

110

111

111

111

100
101

• many states vanishes by interference ⇒ negative (complex) sign problem

8

Statistical Error of MCMC Measurements• There is autocorrelation between successive con gurations2σ(1 + 2τint)
0
2
σ =
M

• σ02 : population variance (determined by the ensemble)• M : number of Monte Carlo steps• τint : autocorrelation time (determined by the MC dynamics)• effective number of independent samples → M/(1 + 2τint)• For systems with negative sign problem

σ02(1 + 2τint)
σ =
M s2
2

• s : average sign (exponentially small for lower temperature, longer time,filarger system)9

Many-body wave function and tensor• Wave function of N-qubit (spin-1/2) system∑| Ψ⟩ =

Cσ1,σ2,⋯,σN | σ1σ2⋯σN⟩

σ1,σ2,⋯,σNN

states → 2N coef cients(Cσ ,σ ,⋯,σ ) should be
1 2
NN
speci ed → memory cost 〜 2• linear combination of 2

• C can be regarded as N-leg (rank-N) tensor
・
・
・

C

σ1

σ2

σ3

σN

• Tensor = multi-dim array = generalization of vectors/matrices• 0-leg tensor → scalar• 1-leg tensor → vector• 2-leg tensor → matrix
• ...

fi

fi

• N-leg tensor → memory/computational cost 〜 exp(N )Tensor representation
NATURE COMPUTATIONAL SCIE
• Tensor is ubiquitous• probability distribution function b

0.8• multi-dim
0.2data0.5

t /T0

P(s1, s2, ⋯, sN )

0

0.8

0.5

• grid data, images

g(x, y) = g(x1, x2, ⋯, xN, y1, y2, ⋯, yN )DNS

y/Lbox

.5

10

0.7
0.5

• (x1, x2, ⋯, xN ) and (y1, y2, ⋯, yN ) are binary number rep. of x and y• e.g) 256x256 image → 256x256 matrix or 2１６
0.3(16-leg) tensor

0.7

χ = 1180.5
0.3• neural network• weight matrix → tensor with many legs

χ = 74

0.7
0.5

・
・
・

・
・
・

W

1.0

1.5

11

Tensor decomposition

・・・• Decomposition of a many-leg tensor into product of small tensors

・
・
・

・
・
・

C

σ1

σ2

σ3

σN

σ1

σ2

σ3

σN

• taking summation of common indices (contraction)• tensor ⇒ tensor decomposition ⇒ tensor network• other representations: Boltzmann machine, path-integral, DNN, etc• Advantage of tensor network• data compression O(exp(N )) → O(N )• freedom in contraction order• highly accurate approximation based on SVD• Tensor network representation• approximation of original tensor (e.g. variational w.f.)• exact representation• quantum states: GHZ state, AKLT state, etc• classical/quantum statistical model, quantum circuit, etc

12Tensor network representation
of
1
国
1国
肆
quantum circuit 10 10 肆
。 >

一

>

。 >

一

𦥯
𦥯

e
e

-

>

1
1 a quantum• Tensor network representationof国
circuit
国
0 >

一

0 >

• all bond dimensions are 2

| 0⟩

⑧

H

-

⑧

一度
。 一度
。
contraction

→

←
→ ←
一

| 0⟩
| 0⟩

-

一

H

H

⑧>
10

>
10④
⑧
国

国

•

一
④

は>

④
は>

1.net
1.net
11
10
一

④

一
一
-

-

一

•

-

朝
朝
11 >

>

10 >
• equivalent to Schrödinger (state-vector)
simulation
學
国
10
11 〉
国學

11 〉

10 >

>

-

• Taking contraction from the initial state (left to right)
-

-

一

>

-

一

一

-

一

• Evaluating an amplitude

| 1⟩

H

| 1⟩ → ⑥

| 0⟩
| 0⟩

H

H

| 1⟩

⑧

⑧

⑥

⑧④

⑧

-

①

-

-

④

①

④

-

①④

一

①

contraction

-

一
-

-

• order of contraction does not change the

-

→
○
○

① は 1 1 は >は 1 1 は >①nal result• freedom in contraction order → possibility to reduce the cost

fi

| 0⟩

13

Optimization of Contraction Order• General guideline for better contraction order• avoid tensors with a large number of legs in the middle or at the end ofcomputation• It is known that nding the optimal contraction order is (at least)#P-hard problem• can'tnd the best solution in a realistic time

• many heuristics have been proposed, c.f.,• Schutski, R., Khakhulin, T., Oseledets, I., & Kolmakov, D., Simple
heuristics for ef cient parallel tensor contraction and quantum circuit
simulation. Physical Review A, 102(6), 1–11 (2020). https://doi.org/
10.1103/PhysRevA.102.062614• Gray, J., & Kourtis, S., Hyper-optimized tensor network contraction.

fi

fi

fi

Quantum, 5, 1–22 (2021). https://doi.org/10.22331/Q-2021-03-15-410

Slicing tensor network• Slicing• for a subset of bonds in the tensor network•

x the bonds to some value• these bonds disappear from the tensor network• for each

xed value• perform contraction independently• take summation on

xed-value patterns at the outermost• Advantage• memory cost of contraction of sliced tensor network becomes smallerfi• contraction of each sliced network can be performed in parallel

fi

fi

14

15

State-of-the-art tensor network simulation• Y. A. Liu et al., Closing the “quantum supremacy” gap: Achieving realTime simulation of a random quantum circuit using a new sunway
supercomputer. International Conference for High Performance
Computing, Networking, Storage and Analysis, SC (2021)• Gordon Bell Prize Winner in 2021

https://awards.acm.org/bell

16

Tensor network methods in statistical physics• “Renormalization" (or “Lagrangian") approach• coarse graining of tensor network representation of the partition function• transfer matrix, tensor network renormalization (TRG), higher-order tensor
network renormalization (HOTRG), etc• “Variational" (or “Hamiltonian") approach• tensor-network approximation of strongly correlated many-body quantum
states• DMRG, PEPS, MERA, etc• “Exact" contraction of tensor network can not be done in two andhigher dimensions• low-rank approximation based on eigenvalue/singular value decomposition• accuracy of approximation is controlled by “bond dimension": D ( or χ)(transmission)

17

Tensor renormalization group• Low-rank approximation based on SVD
computational cost: O(D 5)memory cost: O(D 3)• Improvement of accuracy by considering the “environment" effects
or removing local correlations• second order renormalization (SRG), mean- eld SRG, etc• tensor network renormalization (TNR), loop TNR, Gilt, etc• computational cost increases signi cantlyfifi

Levin & Nave (2007)

18

More advanced tensor-network methods• Improving accuracy• including environment effect• SRG (2009), CTMRG (1996, 2009)• removing local correlations• TNR (2015), loop TNR (2017), Gilt
(2018)• computational cost increases signi cantly• improve accuracy without increasing
complexity• BTRG (2022)• Reducing complexity• TRG using few-leg tensors

comp
Th
gaug
row
the g
flipp
comp
tenso
lows
depi
into
trun
two
chan
in F
resen
a un
Fig.6
glers
gard
comp
twee

• ATRG (2020), CATN (2020)
FIG. 6. (a) Details of the projective truncation made at thefirst step of the TNR iteration; here
two copies ofa projecG. Evenbly(2015)
tor Pu , which is composed of a product isometric and unitary• HOTRG (2012), ATRG (2020), CATN (2020)tensors, are applied to a 2 × 2 block of A tensors. (b) Definition of four index tensor B. (c) Projector Pu is formed fromisometries vL , vR and disentangler u [and their conjugates].
(d) Details of the projective truncation made at the second• Grassmann tensor network (2010, 2021)
step of the TNR iteration. (e) Definition of matrix D. (f)Details of the projective truncation made at the third step of
the TNR iteration. (g) Definition of new four-index tensor
A! , copies of which comprise the coarse-grained square-latticetensor network. (h) Delineation of the different dimensions• Generalization to higher-dimensionsfi• Fermions

By
u† u
cept
rang
cuss
key
mali
indic
to a

v

19

Higher-order tensor renormalization group• Coarse graining in each direction using HOSVD instead of SVD• Squeeze tensors using “isometries"⌘<latexit sha1_base64="6fFfiOA5K65pfNgux7h1VXJdczM=">AAACaXichVFNLwNBGH66vqq+Wi4NF7GpODWzJSEO0sTFkVZLQtPsrsGwX/ajSTX+gJOb4EQiIn6Giz/g4CeII4mLg7fbTYQG72Rmnnnmfd55ZkZzDOH5jD3FpI7Oru6eeG+ir39gcCiZGi57duDqvKTbhu2ua6rHDWHxki98g687LldNzeBr2v5ic3+txl1P2NaqX3d4xVR3LLEtdNUnqrzJDwJRqyZllmVhjLcDJQIyoli2kzfYxBZs6AhggsOCT9iACo/aBhQwOMRV0CDOJSTCfY4jJEgbUBanDJXYfRp3aLURsRatmzW9UK3TKQZ1l5TjyLBHdste2QO7Y8/s49dajbBG00udZq2l5U516DhdfP9XZdLsY/dL9adnH9uYC70K8u6ETPMWektfOzx9Lc4XMo1JdsVeyP8le2L3dAOr9qZfr/DCBRL0AcrP524H5VxWmc7mVmbk/EL0FXGMYQJT9N6zyGMJyyjRuXs4wRnOYy9SSkpLo61UKRZpRvAtJPkT06uMOA==</latexit>• More accurate than TRG in two dimensions• computational cost for 2D: O(D 7)• Works in higher dimensions• computational cost in d-dimensions: O(D 4d−1)
Xie et al (2012)

20

New tensor network algorithms: ATRG and BTRG
• Anisotropic Tensor Renormalization Group (ATRG)• D. Adachi, T. Okubo, S. Todo, Phys. Rev. B 102 054432 (2020)• Works in any dimensions with smaller cost: O(D 2d+1) ≪ O(D 4d−1)• same cost as TRG in two dimensions• More accurate than TRG• Bond-weighted Tensor Renormalization Group (BTRG)• D. Adachi, T. Okubo, S. Todo, Phys. Rev. B 105, L060402 (2022)• Works in two dimensions with the same cost as TRG• More accurate than TRG and HOTRG

21

Central idea of ATRG
• Prevent contracting large tensors• in HOTRG⌘

O(D7 )⌘O(D5 )<latexit sha1_base64="yhmNrzmP8XtfxsaawIJWCG1W4Dc=">AAAB7XicbVDLSgNBEOyNrxhfUY9eBoMQL2FXxeQY0IM3I5gHJGuYncwmY2ZnlplZISz5By8eFPHq/3jzb5w8DppY0FBUddPdFcScaeO6305mZXVtfSO7mdva3tndy+8fNLRMFKF1IrlUrQBrypmgdcMMp61YURwFnDaD4dXEbz5RpZkU92YUUz/CfcFCRrCxUuO2eP1QPu3mC27JnQItE29OCjBHrZv/6vQkSSIqDOFY67bnxsZPsTKMcDrOdRJNY0yGuE/blgocUe2n02vH6MQqPRRKZUsYNFV/T6Q40noUBbYzwmagF72J+J/XTkxY8VMm4sRQQWaLwoQjI9HkddRjihLDR5Zgopi9FZEBVpgYG1DOhuAtvrxMGmcl77J0fndRqFbmcWThCI6hCB6UoQo3UIM6EHiEZ3iFN0c6L8678zFrzTjzmUP4A+fzByf2ji0=</latexit><latexit sha1_base64="6fFfiOA5K65pfNgux7h1VXJdczM=">AAACaXichVFNLwNBGH66vqq+Wi4NF7GpODWzJSEO0sTFkVZLQtPsrsGwX/ajSTX+gJOb4EQiIn6Giz/g4CeII4mLg7fbTYQG72Rmnnnmfd55ZkZzDOH5jD3FpI7Oru6eeG+ir39gcCiZGi57duDqvKTbhu2ua6rHDWHxki98g687LldNzeBr2v5ic3+txl1P2NaqX3d4xVR3LLEtdNUnqrzJDwJRqyZllmVhjLcDJQIyoli2kzfYxBZs6AhggsOCT9iACo/aBhQwOMRV0CDOJSTCfY4jJEgbUBanDJXYfRp3aLURsRatmzW9UK3TKQZ1l5TjyLBHdste2QO7Y8/s49dajbBG00udZq2l5U516DhdfP9XZdLsY/dL9adnH9uYC70K8u6ETPMWektfOzx9Lc4XMo1JdsVeyP8le2L3dAOr9qZfr/DCBRL0AcrP524H5VxWmc7mVmbk/EL0FXGMYQJT9N6zyGMJyyjRuXs4wRnOYy9SSkpLo61UKRZpRvAtJPkT06uMOA==</latexit>• in ATRG<latexit sha1_base64="5LslsaN7jbXUrzYL06zhdVG0JWU=">AAAB7XicbVDJSgNBEK2JW4xb1KOXxiDES5hxzTGgB29GMAskY+jp9CRterqH7h4hDPkHLx4U8er/ePNv7CwHTXxQ8Hiviqp6QcyZNq777WSWlldW17LruY3Nre2d/O5eXctEEVojkkvVDLCmnAlaM8xw2owVxVHAaSMYXI39xhNVmklxb4Yx9SPcEyxkBBsr1W+L1w/nx518wS25E6BF4s1IAWaodvJf7a4kSUSFIRxr3fLc2PgpVoYRTke5dqJpjMkA92jLUoEjqv10cu0IHVmli0KpbAmDJurviRRHWg+jwHZG2PT1vDcW//NaiQnLfspEnBgqyHRRmHBkJBq/jrpMUWL40BJMFLO3ItLHChNjA8rZELz5lxdJ/aTkXZRO784KlfIsjiwcwCEUwYNLqMANVKEGBB7hGV7hzZHOi/PufExbM85sZh/+wPn8ASTsjis=</latexit><latexit sha1_base64="6fFfiOA5K65pfNgux7h1VXJdczM=">AAACaXichVFNLwNBGH66vqq+Wi4NF7GpODWzJSEO0sTFkVZLQtPsrsGwX/ajSTX+gJOb4EQiIn6Giz/g4CeII4mLg7fbTYQG72Rmnnnmfd55ZkZzDOH5jD3FpI7Oru6eeG+ir39gcCiZGi57duDqvKTbhu2ua6rHDWHxki98g687LldNzeBr2v5ic3+txl1P2NaqX3d4xVR3LLEtdNUnqrzJDwJRqyZllmVhjLcDJQIyoli2kzfYxBZs6AhggsOCT9iACo/aBhQwOMRV0CDOJSTCfY4jJEgbUBanDJXYfRp3aLURsRatmzW9UK3TKQZ1l5TjyLBHdste2QO7Y8/s49dajbBG00udZq2l5U516DhdfP9XZdLsY/dL9adnH9uYC70K8u6ETPMWektfOzx9Lc4XMo1JdsVeyP8le2L3dAOr9qZfr/DCBRL0AcrP524H5VxWmc7mVmbk/EL0FXGMYQJT9N6zyGMJyyjRuXs4wRnOYy9SSkpLo61UKRZpRvAtJPkT06uMOA==</latexit>• decompose the local tensor into small pieces and interchange theposition before coarse graining

D. Adachi, T. Okubo, S. Todo, Phys. Rev. B 102 (2020) 054432

T

x0

A

A

α B

m
A
A
E
B 102,
054432(2020)four tensors, A, thThen,
we defineBy1 xy10 PHYSICAL REVIEW
(c)
X
X
T B
e
22
T
C
C
x
A
A
x
0
2
(3)
y1 x1 α = S{T }αα V{T }y1 x1 α ,y1
Ay0 x0 α = U{T }yfr
0xThen, we definey0YYTα four tensors, A, B, C, and D, as
β
th
D
D
T
x2 X(4)TX
x0 ACy1 x2 β = S{Tx1A}ββ U{T }y1 x2 β ,YYn
By1 x1 α = S{T }ααV{T
x
β D
3
=
U
,
(2)
A
y
x
α
{T
}y
x
α
0
0
0
0
D
y2 y1
αT
m
D
D
Y(5)
Y
TXV{T }y2 x3 β ,x2 X Dy2 x3 β =Fy 2 x3Cy1 x2 β = S{T }ββ Uid{T
By1 xx1 α1PHYSICAL= S{T }ααV{TREVIEW(3)}y1 x1 α ,ADACHI,OKUBO,
AND
TODO
B
102,
054432(2020)
(d)
y
(b)
T
β
1
A
[step ADACHI,(a) in Fig.OKUBO,
that the accuracy of the final DfreeDY1]. Note
T
T
ANDYTODOPHYSICALREVIEWlo(e)(a)= V{T }y2 x
Denergy may depend
on
how
to
separatethesingular
value
y
y
x
β
x
C
=
S
U
,
(4)
2
3
0
x1
3
y
x
β
{T
}ββ
{T
}y
x
β
y
1 2
2 1 2
y
0
th
Then,
we
definefour
tensors,
A,
B,
C,
and
D,
as
y
B
matrix in the above
we observe that byAincluding A
x0 G
Gβ D0 step. Indeed,
A
D
(d)
x
TA
Then,
we
definefour
tensors,
A,
B,
C,
A
A
y
T
x0 A
1]. Note that the a0 final free
A the error of the(a) E
the singular 0matrix
B and C,Dy2 x3 β = V{T }y2 x3 β , [step (a) iny1Fig.
x3S inH (5) (2)
H to sep
T
y
T
x
2
(c)
A
A
0
X
energymay
dependon how
[s
=
U
,
A
y0
X
{T }y0 x0 α
energy is minimized.
Such a construction gives us√better y0 x0 α
Cα Bα (d)x=U,Amatrixabove
Indeed,y10 xstep.
{T }y0 x0 αweao
x0 A
A
0α
x2 theMoffree energy thandecomposition,(a) in αFig.
Note
that
the
of in
thethefree
X equal Bweight
X [stepM1].S,AyfinalA accuracy
ADACHI,OKUBO,
AND
TODO
2
E
x
x
the
singularmatrixSinB
X
1 into A and B (orBy1how=toS{Tseparate(3) and C, th
2 depend
the singular ymatrix
C and
If we ondoenergymayD).
x1 αX}αα V{T }y1thex1 α , singular value(f)
D
1 y
α B Y
(c)
Y
1
X
By1 x1 α = SSuch
X Indeed, we observe
energycons{T }αα Va
{T }y
B T
CA Y
C Y(a),matrix
1 x1 αT
not introduce
in A
step
our splittingof singularin thestep.
that isby minimized.
including
x2 truncations
y1above
Then
y
0
T
free
energythanthe(4)
equal weightE
YgivesS C
Y
T free= SC,
U{Terror,xT1 step (b) andD
matrix S, Ttogether with(c) below,
the singulariny1 xB2 βD
andthe{T }ββ}y1 x2 βofD
y1 x1DT matrixTan
x0 theA final
AS{T }ββ
(c)
β
F
C
=
x
β
X
the
singularmatrix
into
AU
and
y
x
β
{T }yB1X
1 2
1 x2 β(identical singularDvalue
directis partial
SVD ofSuch a construction gives us√better
energyminimized.
D Dto theC
D spectrum
x2 C
Y
Y
a
(a)
β
truncations
α, introduce
DyD
T T = ABCD. Inythis
S, of (5) in step (a),freesplittingenergygives
thanthebestequal
D the
β = V{Tdecomposition,}y2 xnotx3 the presentx
2 x3weight
3β
3 sense,
2
(a)
x
y
X
X
y
2
0
Vstep
, to
matrixS,
togethera(e)Dβ Dy2 xforβ ={T }y2 x(b)
x1.
3with
3β
local approximation,
and we think itthe
remains
optimaleven
if A and D
singularmatrixinto
B step(or C
and2 D).
Ifdirection
weDdo
3 D
y
FIG.
Renormalization
of
ATRG
in
y
the
twoy
2
0
(d)
y0
F
x0 A
Y
Y
A
[step
(a)
in
Fig.
1].
Note
that
the
accuracy
of
the
final
free
value spectrum tothere exist truncations in step (a). not introduce
truncations
inlattice
step (a),ouridentical
splittingof singulary1 singulardimensionalsquaremodel.
x xG
G
A
(d)
x0 A
A
Y
Y thetheaccura
T
T
[step
(a)
in
Fig.
1].
Note
that
y20 3
energymay
dependonhowtoseparatethesingular
value
(b)
T
T
=
ABCD.
In
this
sense,
prese
matrix
togetherwithNext, by using
we swap
the S,
bondof B andC step (b) and (c) below, gives an
Dpartial SVD,
D
α
y1 B
x
H
energymay
dependonhowtoseparate
1HBF(e)Indeed,matrixinx thestep.
we observethat byincluding
local approximation,
identical
singularvalue
to thedirectpartialSVDof and we think it[step (b)(c) in Fig.
swap
the
bond
A 1]. In orderA toyaboveα and
0 of Bspectrum
β aboveand1 thevertical(y)directions
alternately
as
in
the
same
way
B
B
E
matrixinthestep.
we(a).
obser
xD
Din stepthesingularmatrix
S
in
B
and
C,
the
error
of
the
final
free
x
1
there exist
truncations
A
T = ABCD.(a)In
thisGsense,Athe present
splittinggives
the bestIndeed,
we define tensor
MTasand x2 bond of C,
yy12 (b)1
x
G
M
A
A
0
(c)
E
as
HOTRG.
The
outlineofone
ATRG
renormalization
step
x
X
the
singularmatrix
S
in
B
and
C,
the
x
1
X(e) !
energyis minimized.
Suchus√partial
better3 eveny2 gives
Next,
by using
SVD, weerr
sw
y1 y0
x2 C
local
and wethinka itconstruction
remains
optimalifyapproximation,(f) Cy0 alongX
(c)
1
X
H
y
direction
is
shown
in
Fig.
1.
H
energyisminimized.
Suchaconstructi
C
C
S,
of
free
energy
than
the
equal
weightdecomposition,
x2 x[step (b) and (c)in Fig. 1]. In orderCy1 xx2 βexist, truncationsA(6)
in step (a).
(d)
GMαβx1 x2 = G By1 x1 αthereAA
0 A First, four-leg
[step
(a
0
tensorT
withbonddimensionχis define
approxiβ
freeenergy
than
the
equal
weightdecoxthesingularmatrix
into
A
and
B
(or
C
and
D).
If
we
do
bond
of
C,
we
tensor
M
and
x
T
1
T
y1
D
D
Next,
by
using
partial
SVD,
we
swap
the
bond
of
B
and
C
M
2
D
M
y
2 a rank-χ matrix by using the partial SVD asenergyβ Dy1 H Y D H Y [stepnot(b)
mated
by
the
singularmatrix
into
A
and
B
(or
C
introduce
truncations
in
step
(a),
our
splittingof
singularα B
and (c)Y in Fig.
to swap the x1 bond of Bx3!
B 1].(f)InY ordermatrix
y2
A
A
and,
by
partial
SVD
of
M
and
truncatingthesingularvaluesnotintroduce
truncations
in
step
(a),
χ
x
matrixS,
together
with
step
(b)
and
(c)
below,
gives
an
M
=
Bour
D
xy32
bond
of
C,
we
definetensor
M
as
and
x
αβx
x
y1 x1
1
D
!
1
2
M
2
M
E
x1
the
sing
y2define new X andb(Fyto χ , weY
asmatrixS,
together
with
step
(b)
and
1
identical
singularvalue
spectrumtothedirectpartial
SVD
of
T
D
T
y
≈
S
U
V
.
(1)
T
D
D
D
1
y
y
x
x
{T
}αα
{T
}y
x
α
{T
}y
x
α
(c)
0 1 0 !
1
0 0X
1 1
X
(f)
alization step of ATRG in y direction for the twoenergycF thisC
identical
value
spectrumtothe di(b)
TxT2 =CABCD.MIn
thesplitting
the
best
=
Bα=1
C
, singulargives(6)
αβx1 x2 sense,
y1 xpresentχ
1 α y1 x2 β
lattice model.
free
en
!
by partial
and trunc(e)T(b)
T T it=and,ABCD.
In thisSVD
sense,
spu
local
approximation,and ywethinkremains
optimalevenofthe
ifM present
DD
1 the truncated SVD) is a method to obtain0
Mαβx1yT
≈
S
U
V
,
(7)
β
The
partial
SVD
(or
sinA
x2
{M}y1 y1 {M}αx2 y1 {M}βx1 y1
DD in step (a).
(e)χ , we define
new
and
Y the
as
local toapproximation,
and X
we
think
it remaexist
y0intruncations
Y isometries
Yby
xFIG.
G stepthere
GRenormalization
A
A
1.
of
ATRG
y
direction
for
the
two0
y1
5
largestsingularvaluesandcorresponding
using
not
intr
y) directions
alternately as in the same waycost: O(D )x3byG
and, bypartialSVD
MG
andtruncating
thethesingularvaluesthere
exist
truncations
step• computational
usingofpartial
SVD,
we swapbond of
BinandC (a).
Adimensional
A " latticeχ
yNext,
0 Arnoldi
squaremodel.
2 xtheymethodorthelow-rank
approximation
based
on matrix
a
!
1
H
H
D
Xαx2 y1 = S{M}y1 y1 Uto(8)outline of one ATRG renormalization stepD{M}αxy1 , define
χ
,
new
X
and
Y
as
Next,
by
using
partial
SVD,
2 we
[step
(b)
and
(c)
in
Fig.
1].
In
order
to
swap
the
x
bond
of
B Swe swap
1F
3
y
M
≈
U t
FIG.
1.
Renormalization
step
of
ATRG
in
y
direction
forthetworandomized
technique[32,33].
By
using
the
partial
SVD,
the
1
αβx
x
{M}y
H
identica
1 2
1 y1 {MHs shown in•Fig.
1.
O(D)
memory:x[steporder to sw
define tensor
M(b)as and (c) in Fig. 1]. In
and x2 bond of C, we
1 "M model.
M
y2
χcost
y1
latticecomputation
for calculating
K largest singular valuesof = A
(b)as
!
and
the
vertical(y)directions
alternately
in
the
same
way
Y
=
S
V
.
(9)
ensor T with bond dimensional
dimension squareχ is approxix
βx
y
{M}y
y
{M}βx
y
1
1
1
1
1
1
bond
of
C,
we
definetensorandx" TMT as1
M (f)2
M dimensions
!
y2M
M
≈
S
U
V
,
(7)
(e)
•
generalizationtohigherisstraightforwardanbyNmatrix
(K
"
M,
N)
can
be
reduced
to
O(KMN
),
αβx
x
{M}y
y
{M}αx
y
{M}βx
y
local
ap
1 2
Xαx2 y(6)= S{M}yU
as HOTRG. The outline of one ATRG renormalization
χ matrix by using the partial SVD as1 !
1 y1Myαβx1 x2 =1 1stepBy1 x21 α1 Cy1 x2 βy,10 1
(f)
which
is
significantly
smallerthanO(min(M,
N
)MN
),
the
1
Then,
we
renormalize
thehorizontal
two
bonds
into
one
there
ex
directionis shown
in Fig. 1.
and the vertical(y) directions
alternatelyalong
as inythex0 G Mαβx1 x2 =G "By1 x1 αCy1 x2
A
A "
χ
T
T same way
y1DD!
computation
cost
of
the
standard
full
SVD. It has(8)been shownby usingsqueezersE
and
F
[step
(d)
and
(e)
in
Fig.
1].
We
S{M}yNexV
X=
S
U
First,
four-legtensorT
withbonddimensionχisapproxias
HOTRG.
The
outlineofone
ATRG
renormalization
step
βx1 y1 =αx2 y1 Phys.
y1 , y1(2020) Y
T {M}yRev.
1 y1 F1 y1 {M}αxD.
B 2102054432DD Adachi, T. Okubo, S.TTodo,
S{T }ααU{T }y0 x0 αV{T }y1 x1 α .
(1)
x1 ≈
H
H y1 TRG
that
by
using
this
technique,one
canperformtheoriginal
call
themsqueezers
since
they
are
not
necessarily
isometries,[step(balong y direction is shown in Fig. 1.
mated by a rank-χ matrix by and,usingbythe
partialSVDSVDof"
asM and truncating the singular values
partialpα=15x)computation
cost
[11].
In
the
present
calculation,
with
O(χ
Then,
we
renormalize
thehorizounlikeBy applying
squeezer
E
(F
)
to
b
and
x
1
and,
by
partial
SVD
of
M
and
truncating
M
Y
=
S
V
.
(9)
First, four-leg tensor T withbondconventional
dimension χHOTRG.
is approxi2
to
χ
,
we
define
new
X
and
Y
as
M
βx1 y1{M}y1 y1 {M}βx1 y1 y2χ!
5FIG. 1. Renormalization step of ATRG in y direction for the two-

α OKUBO, AND TODOADACHI,x2 X
X

Ay0 x0 α = U{T }y0 x0 α ,(2)B

One Renormalization Step of ATRG

23

Benchmark Test of ATRG in 2D• Free energy of square lattice Ising modelcomparison at same bond dimension
10-5

TRG
ATRG
HOTRG

comparison at same computational cost

10

TRG
ATRG
HOTRG

-5

10-6

δf

δf

10-6
10-7

10-7
10-8

10-8
10-9
10

20

50

100

106

107

108

109

τ

χ
D
<latexit sha1_base64="gtucZmzsqVUGxXRQ1mOyFbKRsEs=">AAACZHichVFNSwJBGH7cvswsLQmCICQxOsloQdJJqENHP/IDTGR3G2tx3V12V8GkP1DXokOngojoZ3TpD3TwDwTR0aBLh17XhSip3mFmnnnmfd55ZkYyVMWyGet6hJHRsfEJ76Rvyj89EwjOzhUsvWnKPC/rqm6WJNHiqqLxvK3YKi8ZJhcbksqLUn2rv19scdNSdG3Xbhu80hAPNKWmyKJNVGa7GoywGHMiPAziLojAjbQevMUe9qFDRhMNcGiwCasQYVErIw4Gg7gKOsSZhBRnn+MYPtI2KYtThkhsncYDWpVdVqN1v6blqGU6RaVukjKMKHtid6zHHtk9e2Efv9bqODX6Xto0SwMtN6qBk4Xc+7+qBs02Dr9Uf3q2UUPS8aqQd8Nh+reQB/rW0UUvt5mNdlbYNXsl/1esyx7oBlrrTb7J8OwlfPQB8Z/PPQwKiVh8LZbIrEdSSfcrvFjEMlbpvTeQwg7SyNO5HKc4w7nnWfALIWF+kCp4XE0I30JY+gSMmIm+</latexit>• ATRG outperforms TRG at the same bond dimension⌧⌘(

1010 1011 1012

D5
D7

TRG and ATRG
HOTRG

<latexit sha1_base64="pJxdtBZkOPLCnS7faCubKAJn2NA=">AAACuXichVFNT9RAGH6oH0D9YNWLiZfGDcTTZhYkGEwMRhO5CcsukFDcTLsvy0h3WtrpRmz2D/AHPHjSxBjjP/DqxT/ggZtX4IiJFw++7TYaJerbdOaZ532fd56Z8aJAJUaIgxHrzNlz50fHxu0LFy9dnqhcubqahGnsU8sPgzBe92RCgdLUMsoEtB7FJHteQGvezoM8v9anOFGhbpq9iDZ7sqvVlvKlYapdabpGpo5Lu6nqO65HXaUzn/slA/vhk1lnynENPTNZs/HIkbrj3GcwcF3Ozf3MLT7OSdsl3Sml7UpV1EQRzmlQL0EVZSyFlbdw0UEIHyl6IGgYxgEkEv42UIdAxNwmMuZiRqrIEwawWZtyFXGFZHaHxy6vNkpW8zrvmRRqn3cJ+I9Z6WBSfBbvxIn4JN6LQ/H9r72yokfuZY9nb6ilqD2xf33l239VPZ4Ntn+p/unZYAt3Cq+KvUcFk5/CH+r7z1+crMw3JrMp8Vocs/9X4kB85BPo/lf/zTI1XsLmB6j/ed2nwep0rT5Tm16+XV24Vz7FGG7gJm7xfc9hAYtYQov3/YAvOMSRddeS1rb1dFhqjZSaa/gtrOQHbLWo9A==</latexit>• ATRG outperforms HOTRG at the same computational cost

D. Adachi, T. Okubo, S. Todo, Phys. Rev. B 102 (2020) 054432

24

Contracting arbitrary tensor networks• Combination of• contraction of tensor network with optimized order• low-rank approximation using SVD (if bond dimension exceeds D)• avoid many-leg tensors and keep MPS form of three-leg tensors

2

FIG. 1. Detailed process of contracting a tensor network with 5 nodes, each of which is a four-way tensor, as sketched in Fig. 2 of maintext. The scissor symbols in the figures indicate applying SVDs on the matrices unfolded from the tensors. The tensors in the step (1) are theF. forming
Pan, P.
Zhou,
S. pentagon.
Li, P. Zhang,
Phys.
Rev.
125,transformed
60503 (2020)
original four-way tensors connected to each other,a fully
connected
The step (2)shows theMPS Lett.
representationfrom (1); the arrow indicates contracting two MPSes, that is, annihilating one color. As an example, the green MPS and the purple MPS iscontracted to a longer purple MPS in (3). The steps (3) (5) show the swap operation between two tensors in the purple MPS. To accomplish

25

Markov-chain Monte Carlo• Various update algorithms• local, cluster, worm, event-chain, hybrid (HMC), etc• extended ensemble methods• Pros• consistent in long-time limit (if balance condition and ergodicity aresatis ed)• computational cost increases only linearly (in many cases)• trivial parallelization• Cons• slow convergence of statistical error ( ∼ 1/ M )• effective number of samples decreases as auto-correlation increases2σ(1 + 2τint)
0
σ2 =
M

• negative sign for frustrated quantum magnets, fermions, real-time

fidynamics

26

Tensor network renormalization group• Many variants• tensor renormalization group (TRG) by Levin and Nave (2007)• higher-order tensor renormalization group (HOTRG) by Xie et al (2012)• tensor network renormalization (TNR) by Evenbly and Vidal (2015)• anisotropic tensor renormalization group (ATRG) by Adachi et al (2020)• bond-weighted tensor renormalization group (BTRG) by Adachi et al(2022)• Pros• (super-?)exponentially accurate for large bond dimension D• Cons• systematic error (bias) due tonite bond dimension D• computational cost increases rapidly as Dαfi• few methods for higher dimensions

27

Approximate contraction of tensor network• “Exact" contraction of tensor network can not be done in two andhigher dimensions• low-rank approximation based on singular value decomposition• systematic bias due to low-rank approximation• accuracy of approximation is controlled by cutoff (or bond dimension: D)• convergence is fast but not smooth (even not monotonic) in manycases

• How can we eliminate the systematic bias in approximate
contraction?
• combination with Markov chain Monte Carlo?

28

Advances in Markov chain Monte Carlo• Representation (de nition of “con gurations” and “weighs”)• path integral representation for quantum Monte Carlo (1976), Bayesianinference (1990), tensor network representation• Choice of ensemble• extended ensemble method: multicanonical MC (1991, 2001),
exchange MC (1996), lifting (2000)...

• Generation of set of candidate con gurations• non-local (cluster) updates: Swendsen-Wang (1987),
Hamiltonian MC (1987), loop (1993), worm (1998)...

• Choice of transition kernel (probabilities)• Metropolis, heat bath (Gibbs sampler), over-relaxation (1987),
irreversible kernel (2010), event-chain (2013)...

• Algorithm for generating a con guration according to transition
probabilities

fi

fi

fi

fi

• N-fold way (rejection free) (1975), Walker’s method (1977, 2019),
order-N algorithm (1995, 2009)...

29

Agenda• Introduction• Markov-chain Monte Carlo and negative sign problem• Tensor network representation and approximate contraction• Sampling approach in tensor network contraction• Projector formulation of tensor network method• Sampling projectors• Monte Carlo in tensor network representation• Sequential Monte Carlo in tensor network representation• Markov-chain Monte Carlo in tensor network representation• Results• SummaryProjector formulation of tensor network
methods• HOTRG (Xie et al 2012)30

Projector formulation of tensor network
methods

31• TRG (Levin-Nave 2007)• original truncation based on SVD

based
rmca.fitrmca.fi
on

on

.

旦⑤旦

STD
based

on

STD

on

旦

旦

旦⑤旦 @
→

Project r
。

や

?
旦

三

旦
=

や

D

旦

R

R

R

旦

-

@

R

R
-

R
i

→

?

R

-

-

R
-

旦

%

?

一

i

-

%

→

一

• truncation based on projector

Project r
。

旦

→

→

婿

題旦

?

P

題旦

R

PR

18 も
旦

→

=

R

D

きが
きが𣳾
三

-

で

R
一

R

R
-

=

D

き

で

R

旦

R

一

18 も

が

が

- . -。

STD

が

-

_

=

き

旦

旦
③旦

→

32

Levin-Nave TRG in projector formulation• 4 × 4 square lattice case (N = 16)

-

0

-

-

0

8

-

-

8

L

|L
|

-

-

1

-

も𤇾
も𤇾
1

!

羅籩

憗
燻
礴※
唯
嗹
第
※
第店 明※
夔
簿
港
※
※
夔:
店
港
※簿鵼毱
※

爕
:
癭
爕
:
癭
!

!
燻怔!が
礴越滅憗
趙

趙

の
の
Y

o
o

Y

0- _ -☆
☆
.

0- _ -

が州兼 W W † 牧華
州兼θ θ 牧華.

• approximate partition function• contraction of tensor network of 2N initial tensors and (N − 4)projectors• contraction graph of depth ∼ log N

Projector formulation of tensor network
methods

33• ATRG (Adachi et al 2020) and CATN (Pan et al 2020)
PHYSICAL REVIEW B 102, 054432 (2020)

ADACHI, OKUBO, AND TODO

T

T

x0

A

A

X

X

Y

Y

x1
β D
y 2 x3

D

α

x2
T

y1

T

(a)
y0
x0 A

A

α B

B

x1
C

C

β D

D

y1

x2

y2

(d)
A

(c)

E

Y

Y

D

D

F

(b)
A

y0
x0 G
y1 H

A

M

M

D

D

A
X

X

x3

(e)
G
H

x1

y2

(f)
T

T

FIG. 1. Renormalization step of ATRG in y direction for the twodimensional square lattice model.

and the vertical (y) directions alternately as in the same way
as HOTRG. The outline of one ATRG renormalization stepalong y direction is shown in Fig. 1.
First, four-leg tensor T with bond dimension χ is approximated by a rank-χ matrix by using the partial SVD as
Ty0 y1 x0 x1 ≈

χ
!
α=1

2

Then, we define four tensors, A, B, C, and D, asy0

S{T }ααU{T }y0 x0 αV{T }y1 x1 α .

(1)The partial SVD (or the truncated SVD) is a method to obtainlargest singular values and corresponding isometries by using
the Arnoldi method or the low-rank approximation based on a

Ay0 x0 α = U{T }y0 x0 α ,(2)

By1 x1 α = S{T }ααV{T }y1 x1 α ,(3)Cy1 x2 β = S{T }ββ U{T }y1 x2 β ,(4)

Dy2 x3 β = V{T }y2 x3 β ,(5)[step (a) in Fig. 1]. Note that the accuracy of the final freeenergy may depend on how to separate the singular value
matrix in the above step. Indeed, we observe that by including
the singular matrix S in B and C, the error of the final freeenergy is minimized. Such a construction gives us√better
free energy than the equal weight decomposition, S, ofthe singular matrix into A and B (or C and D). If we do
not introduce truncations in step (a), our splitting of singularmatrix S, together with step (b) and (c) below, gives anidentical singular value spectrum to the direct partial SVD of
T T = ABCD. In this sense, the present splitting gives the best
local approximation, and we think it remains optimal even if
there exist truncations in step (a).
Next, by using partial SVD, we swap the bond of B and C[step (b) and (c) in Fig. 1]. In order to swap the x1 bond of B
and x2 bond of C, we define tensor M as
FIG. 1. Detailed processof contracting a tensor network with 5 nodes, each of which is a four-way tensor, as sketched in Fig. 2 of main!
By1in
Cy1 xfigures,(6)
αβx1 x2 =x1 αthetext. The M
scissorsymbols
indicate applying
SVDs on the matrices unfolded from the tensors. The tensors in the step (1) are the
2βy1original four-way tensors
connected to each other, forming a fully connected pentagon. The step (2) shows the MPS representation transformed(1);
the ofarrow
indicates
contracting
two values
MPSes, that is, annihilating one color. As an example, the green MPS and the purple MPS is
and, by frompartialSVD
M andtruncating
the singularlongerto χ , wecontracted
define newtoXa andY aspurple MPS in (3). The steps (3) (5) show the swap operation between two tensors in the purple MPS. To accomplish
swapping, we contract two purple tensors first, then apply the SVD on the contracted tensor as shown in (3) and (4). Note that in (5), weχkeep the canonicalform of the MPS. In steps (6) (8), the swap operation is repeated until two tensors connecting the same pair of colors are!
Mαβx1 x2to≈adjacent
S{M}ypositions.
U
V steps,(7)
switched
1 y1 {M}αx2 y1In{M}βx1 y1(9) and (10), we finish the merge step by contracting two tensors indicated by the arrow, producing
y1
a thick bond between
redandpurpleMPSes.
The
steps (10) (14) represent the merge between the purple-and-brown MPS pair and the"
purple-and-blue
pair.
depict the procedure mentioned above repeatedly until a scalar left in the end of the wholeXαx2 y1 =MPS
S{M}yUThe steps, (14) (22) (8)1 y1 {M}αx2 y1
contraction process.
"Yβx1 y1 = S{M}y1 y1 V{M}βx1 y1 .
(9)Then, we renormalize the horizontal two bonds into one
process andis illustrated
tensordiagramby usingThesqueezersEF [step (d)using
and (e)in Fig.
1]. We notations in Fig. 2. To ensure that summing over the index i results to anotherMPS,
we first
dotheyswap
to switchisometries,
the indices i and kcall themsqueezers
since
areon
notA
necessarily
unlike conventional HOTRG. By applying squeezer E (F ) to
X
X
(23)
A and X (Y and D), we obtain new tensor G (H) asai jk =
a(1)
a(1)
a(2)
a(3)i ,
(3)j↵ a↵ik ⇡
j↵ b
↵k b
!
↵
↵,
$ =
$,
G
A
X
E
(10)

Project r
。

旦⑤旦

@旦

旦

→

R

R

R

R

-

→

旦

%

?

i

-

一

-

③旦

→

P
旦
旦
?
題
Projector
formulationof tensor networkきProject旦? 旦 題旦
18 もP
methods
き
旦
三

。r

R

R

D

きが
きが 𣳾

で

-

旦 や

三

=

R

-

や

旦

=

R

D

で

R

R

R

一

=

R

一

=

18 も

• ATRG (Adachi et al 2020) and CATN (Pan et al 2020)• leg swap based on SVD

D

D

婿
婿

が

→

→

が
が

𣳾
が

- . -。

- . -。

STD

-

-

_

_

STD

• leg swap based on projector

iolgri嬴
注が嬴
がiolgri

→

-

-

壇

→

→

→

注

• Any tensor network renormalization methods can be reformulatedusing projectors (?)

壇

34

arXiv:1507.00767

Unbiased Monte Carlo for the age of tensor networks
Andrew J. Ferris

ICFO—Institut de Ciencies Fotoniques, Parc Mediterrani de la Tecnologia, 08860 Barcelona, Spain(Dated: July 6, 2015)
A new unbiased Monte Carlo technique called Tensor Network Monte Carlo (TNMC) is introduced basedon sampling all possible renormalizations (or course-grainings) of tensor networks, in this case matrix-productstates. Tensor networks are a natural language for expressing a wide range of discrete physical and statisticalproblems, such as classical and quantum systems on a lattice at thermal equilibrium. By simultaneously sampling multiple degrees of freedom associated with each bond of the tensor network (and its renormalized form),
we can achieve unprecedented low levels of statistical fluctuations which simultaneously parallel the impressive
accuracy scaling of tensor networks while avoiding completely the variational bias inherent to those techniques,
even with small bond dimensions. The resulting technique is essentially an aggressive multi-sampling techniquethat can account for the great majority of the partition function in a single sample. The method is quite generaland can be combined with a variety of tensor renormalization techniques appropriate to different geometries anddimensionalities.
INTRODUCTION

the mid-1900’s onwards, we have seen numericaldent samples, reducing the denominator. Importa
niques to reduce autocorrelation include cluster up
avoid the well-known slowing down in critical probl

36

Random choice of projectors (Ferris 2015)• Choose D vectors from R (=D 2) basis vectors
R
N
=
p
•
(D) candidate projectors• basis vectors are generated by SVD and chosen randomly according to its
singular value• Stochastic projector is generated as a linear combination of D
projectors• such that on average,

5

1
WθWθ† = E
Np ∑
θ
(a)

(b)

4
(c)

FIG. 1: (a) A two-dimensional partition function expressed as a tensor network. Boundary-MPS contraction proceeds by combining the top
row of tensors with the next (red circles), and then (b) projecting the combined horizontal bonds to a subspace of maximal dimension D(triangular tensors). Typically, the projectors are chosen to maximize the fidelity of resulting effective state of the upper-half of the system,
outlined in blue. (c) After repetition, the end result is a tensor network that is contractible with cost linear in system size.

Graphically, this relation is:
1 X
|S|
s2S

partitionfunction of six-vertex model (16x16)
fixed particle number, D). How to sample efficiently
and directly from this distribution is not immediately obvious, butFIG. 3: Calculations of the partition function of the six-vertex model
it is possible to write the probability distribution itself as aon a 16x16 lattice (where Zexact ⇡ 7.05 ⇥ 1056 ), using standard=(5)
simple matrix-product state56from which perfect sampling canMPS truncation and 104 tensor network Monte Carlo samples with
FIG. 2: One of the approximately
7.05 ⇥ 10 loop configurations
be performed [32, 35]. See the Appendix for details and the37

Summary of arXiv:1507.00767• “Perfect" sampling of tensor network contraction• can be applied to other tensor network renormalization variances (e.g.
HOTRG, ATRG)• only works fornite systems• unbiased? → yes (for partition function)• NB: biased for free energy, expectation value of physical quantities• suffers from exponentially large variance• proposed method = sequential Monte Carlo without resampling• variance of weight of “walkers" is multiplicative and increasesfiexponentially38

Agenda• Introduction• Markov-chain Monte Carlo and negative sign problem• Tensor network representation and approximate contraction• Sampling approach in tensor network contraction• Projector formulation of tensor network method• Sampling projectors• Monte Carlo in tensor network representation• Sequential Monte Carlo in tensor network representation• Markov-chain Monte Carlo in tensor network representation• Results• Summary39

Sequential Monte Carlo• aka• green's function Monte Carlo• transfer matrix Monte Carlo• population Monte Carlo• Monte Carlo• particle

lter

lter

• bootstrap

lter• SISR (sequential importancesampling with resampling)https://www.researchgate.net/ gure/Sequential-Monte-Carlo-scheme_ g2_322302619• Central idea of SMC• approximate probability distribution by (weighted) ensemble of particles• by using Markov chain

fi

fi

fi

fi

fi

• control variance by resampling40

A simple example
x0 = 1
xk = ξk xk−1 (k = 1,2,⋯, n)• ξk is sampled independently randomly from uniform distribution between 0
and 2• expectation value: ⟨xn⟩ = 1• variance increases rapidly for large n

41Resampling• Simple sequential importance sampling becomes unstable for large
steps• weight of each walker is updated randomly by weight factors: W = w1w2w3⋯• random walk diffusion in logarithmic scale• weight degeneracy: weight variance (discrepancy between weights) grows
exponentially and only a few walkers dominate• Resampling is necessary to stabilize the algorithm• resampling:

Pi ≃

∑
k

Wk δi,ik ⇒δ∑ i,ĩkk• after resampling, all walkers share the same
weight:∑k

Wk /Nw

Li-Stattar-Sun (2012)

42

Effect of resamplingx0 = 1
xk = ξk xk−1 (k = 1,2,⋯, n)• ξk is sampled independently randomly from uniform distribution
between 0 and 2
2

10without resamplingwith resampling

1

σ2

<xn>

1.5

without resamplingwith resampling

1

0.1

0.5

0.01

0

0.001

10
n

10
n

43

Tensor network sequential Monte Carlo• Markov chainp(θ1, θ2, ⋯, θn) = p(θ1)p(θ2)⋯p(θn)• sequential selection of projectors• Walker• each Markov chain sequence• a number of Markov chains run simultaneously• Weight• partition function with present projector con guration• Resampling• necessary for controlling the variance between walkers• Systematic error in physical quantities• ∼ 1/M(M: number of walkers)• Memory costfi

• M walkers should be simulated simultaneously

44

Markov-chain Monte Carlo approach• Determine projector candidates from SVD during the conventional(deterministic) TRG• projectors becomes independent with each other and can be sampledindependentlyp(θ1, θ2, ⋯, θn) = p(θ1)p(θ2)⋯p(θn)• (exact) tensor network representation of partition function

Z=

∑
{θi}

g(θ1, θ2, ⋯, θn)p(θ1, θ2, ⋯, θn) =∑{θi}g(θ1, θ2, ⋯θn)p(θ1)p(θ2)⋯p(θn)• Sample projectors {θi} using Markov-chain Monte Carlo• propose new θi according to p(θi)• Metropolis update with P = min(1,g(θ1, θ2, ⋯, θ′i , ⋯, θn)/g(θ1, θ2, ⋯, θi, ⋯, θn))
• update of weights is O(log N ) and includes matmul only
• SVDs are required during initialization stage only

• Physical quantities• can be evaluated by using impurity tensor technique (without systematic
￼

bias)

45

Comparison with Levin-Nave TRG + impurity tensor• Square-lattice Ising model (L = 8)
-0.8

-1.2

specific heat-1

energy

1.2

D=1
D=2
D=3
D=4
exact
TNMC (D=4)

-1.4
-1.6
-1.8
-2

0.6
0.4
0.2

1

1.5

2
T

2.5

0

3

4

1.5

2
T

2.5

2
T

2.5

3

0.9
0.8

magnetization2

3
2.5
2

D=1
D=2
1.5
D=3
D=4
1
exact
TNMC (D=4)
0.5
1
1.5

1
1

3.5

energy2

D=1
D=2
1
D=3
D=4
exact
0.8
TNMC (D=4)

D=1
D=2
D=3
D=4
exact
TNMC (D=4)

0.7
0.6
0.5
0.4
0.3
0.2
0.1

2
T

2.5

3

0

1

1.5

3

46

Comparison with Metropolis algorithm• Square-lattice Ising model (L = 8, T = Tc, Nmcs = 8192)• statistical error is smaller by orders of magnitude• statistical error decreases exponentially as

D is increased

σ02(1 + 2τint)
σ =
M
2

10-2
Metropolis (E)(m2)
TNMC (E)(m2)

statistical error

10-3

10-4

10-5

10-6

4

6

8
D

10

12

47

Ising model in imaginary external

eld
G Mussardo et al

J. Phys. A: Math. Theor. 50 (2017) 484003• Square lattice Ising modelH=−

Yang-Lee zeros on complex plane of fugacity z

σσ −h
σ
∑ i j
∑ i
⟨i,j⟩

i

• pure imaginary externaleldz = e −2βh = − 1

h = iπ/2β• non-positive Boltzmann weight (m: total magnetization)of the fugacity z.

Figure 1. Distribution of the Yang–Lee zeros for the Ising model in the complex plane

W = e β ∑ σiσj × (−1)m/2

accumulate in certain regions or curves of the complex plane, with their positive local density
η(z) which changes by changing the temperature; if at a critical value Tc, the zeros accumulate
and pinch a positive value of the real axis, this is what may mark the onset of a phase transition.
As we are going to discuss extensively through the rest of the paper, the pattern of zeros of
grand-canonical partition functions can be generally quite interesting and this study alone is asource of many stimulating physical and mathematical questions. If the study of the patternsof Yang–Lee zeros is then the first topic of this paper, the Yang–Lee model (and its zeros!) is
our second main topic. In order to introduce such a model and present the work of this paper
1
in its proper perspective, we need to talk about the pattern of zeros of just one particular statistical system: the Ising model. In [3] Yang and Lee showed that for ferromagnetic Ising-likemodels, independently on the dimensionality and regularity of the lattice and also largely
independently on the nature of the couplings, the zeros of the Ising model lie on the unit circle5 in the complex plane of the variable z = e−2βh (where β = 1/(kT) and h is the external
magnetic field): posing z = eiθ , they have the following structure (see figure 1):(m: total magnetization)• Standard Markov chain Monte Carlo suffers from severe negative
sign problem
1
L=4
8
16
32

average sign

0.6

• for T > T0.1
c the zeros are placed along a ‘C’, namely a symmetric arc around θ = π whoseedges are at ±θ0 (T);• at T = Tc these edges move to the real axis and pinch it;• for T < Tc the zeros densely cover the entire circle.

0.2
0

1

1.5

2
T

2.5

Kortman and Griffiths [9] were the first to notice that the density of the Yang–Lee zeros of
the Ising model nearby the edges ±θ0 (T) gives rise to a problem which has its own interest
since such a0.01
density presents an anomalous behaviour with a scaling law ruled by a critical
3 σ
1
1.5
2
2.5
3
exponent
T T> Tc .

η(θ, T) ∼ |θ − θ0 (T)|σ ,
fi

0.4

fi

average sign

0.8

(1)

48

Ising model in imaginary external

eld• Our proposed method also has negative signs for small D• D = 2 results are almost similar to the standard method• NB: negative signs can appear for small D even if the original model is freefrom negative sign

• However, average sign is improved drastically as we increase D

L = 16

L = 32

1

1

0.8

0.8

0.8

0.4

standard
D=2
D=3
D=4
D=6
D=8

0.2
0

1

1.5

0.6
0.4
0.2

2
T

2.5

0
3 1

0.6
0.4
0.2

1.5

2
T

2.5

0
3 1

1.5

fi

0.6

average sign

1

average sign

average sign

L=8

2
T

2.5

3

49

Summary• Markov chain Monte Carlo in tensor network representation• reducing statistical error using approximate tensor network contraction• removing systematic bias by sampling singular vectors (projectors) using
MCMC→ avoid divergence of statistical error, negative signs and systematic bias

• Computational complexity of one Monte Carlo update• O(D α N log N )• matmul only (no SVD) during MCMC sampling → ideal for modern GPGPU or
HPC• Similar formulation is possible for sequential Monte Carlo withresampling• advantage in calculating free energy/partition function• Applications: HOTRG (2012), ATRG (2020), quantum spin models,fermions, lattice QCD, quantum circuits, etc

S. Todo (in preparation)