京都大学

数理解析研究所(RIMS) 研究集会「次世代計算科学の基盤技術とその展開」
RIMS研究集会報告書 No.1848 (2012) 83-92.

多体問題におけるテンソルネットワーク法

Tensor network schemes for many-body problems
京都大学 大学院情報学研究科

原田健自

Kenji Harada
Graduate school of informatics, Kyoto University

1 序
物質の性質（物性）は多数の要素が相互作用している系の巨視的な性質として記述される多体問題として研究されてきた．しかし，多体問題は解析的にも数値的にも難しい問題として知られている．例えば，物性研究における代表的な多体問題として，量子多体系の基底状態計算や古典系の自由エネルギー計算等があるが，その定義通りの計算量は系のサイズに対して指数関数的に増大する．そこで，解析的および数値的研究において，さまざまな近似的手法が提案され適用されてきた．数値的計算手法においても，近年，テンソルネットワークを用いた一群の手法が，従来法の強力な一般化として登場し，注目を集めている．例えば，代表的なテンソルネットワークとして，Multi-scale entanglement renormalization (MERA)[1], Projected entangled

pair state (PEPS)[2] 等が提案されている．本稿では，これらのテンソルネットワークを量子系の基底状態計算や古典系の自由エネルギー計算に導入する手法，更に，その計算に必要なテンソル縮約や Higher-ordersingular value decomposition(HOSVD)[3] について紹介する．

2 テンソルネットワーク
本章では，物性研究で扱われる多体問題（量子系の基底状態，古典系の自由エネルギー）を通して，テンソルネットワークとは何かを解説し，更に，代表的なテンソルネットワーク (MERA, PEPS) の紹介を行う．

2.1 量子系の基底状態
量子力学に従う物質の状態は量子状態と呼ばれヒルベルト空間で記述される．例えば，結晶中では格子点上に原子が規則的に並んでおり，その原子に拘束された電子が結晶の磁気的特性を決めている．電子の磁気的性
質は，S=1/2 スピンと呼ばれる量子的なスピンで記述され，量子スピン状態は次のようなアップとダウンという２つの古典スピン状態の重ね合わせになっている．

|Ψ ⟩ = a| ↑ ⟩ + b| ↓ ⟩(|a|2 + |b|2 = 1).

(1)

つまり，量子 S=1/2 スピンは２次元のヒルベルト空間のベクトルとして記述される．更に，N 電子系の量子スピン状態は各量子スピン状態の直積を用いて以下のように記述される．

|Ψ⟩ =! !

s1 =↑,↓ s2···!
sN

1

Ts1 s2 ···sN |s1 s2 · · · sN ⟩.

(2)

この係数 Ts1 s2 ···sN は確率振幅と呼ばれ，今の場合，2N 個ある．このように，一般的に量子系は重ね合わせの状態を取ることができることから，その状態を記述するためには，システムサイズに対して指数空間的に増大する要素数をもつベクトルが必要となる．
各電子同士は相互作用をしており，例えば，代表的な二体スピン相互作用として反強磁性ハイゼンベルグ相互作用がある．そのハミルトニアン演算子は次のように定義されている：(格子点 i と格子点 j 上のスピンの反強磁性ハイゼンベルグ相互作用)⎡⎤⎛⎞
1 00 0| ↑↑ ⟩⎢⎥⎜⎟⃗i · S⃗j = 1 ⎢0 −2 2 0⎥ ⎜| ↑↓ ⟩⎟ .
Hij = S⎣⎦⎝| ↓↑ ⟩⎠
4 0 2 −2 00 00 1| ↓↓ ⟩

(3)

!

(4)

また，全エネルギーは近接原子上の電子同士の相互作用エネルギーの和で次のように定義される．

H=

Hij .

ij

このように全系のハミルトニアンは局所相互作用の和で定義されることから，一般的に全系のハミルトニアンは疎な演算子であることが特徴である．
ハミルトニアンの固有状態が温度 T で観測される確率は，そのエネルギー固有値 E を用いたボルツマン重み exp(−E/T ) に比例する．従って，低温で重みを増す最低エネルギー固有状態が物性研究では最も関心が高く，基底状態と呼ばれている．全系のハミルトニアンは疎であることから，ランチョス法などを用いて基底状態を求めるということはよくなされている．しかし，システムサイズに対してヒルベルト空間のサイズが指数関数的に増大するため，かなり小さな系の計算（厳密対角化と呼ばれる）しかできない：例えば，39 個のスピン系では，239 ≈ 約６８９億次元の問題を扱う必要がある．

このように量子系の基底状態計算はかなり難しい問題であるが，１次元系に限れば密度行列繰り込み群(Density matrix renormalization group, DMRG) と呼ばれる非常に成功した手法がある．DMRG の特徴は，最低エネルギー状態を探索する空間を行列積状態 (Matrix product state, MPS) と呼ばれる状態に制限した点にある．DMRG では MPS 状態の中で最低エネルギーをもつ状態を計算する．MPS 状態は確率振幅が次のような行列積で定義される量子状態である．

Ts1 ···sN ≡ ⟨s1 · · · sN |Ψ⟩ ≡ Tr[M (s1 ) · · · M (sN )],(5)ここで，M (si ) はアップとダウンの２種類の χ 行 χ 列の行列である．行列積の順番は１次元系の格子点の順番に取る．つまり，隣接格子点同士の行列積を取っている．このような定義はもっと一般化することができる．それを行うために，si （アップとダウン）もインデックスと考え，M (si ) を３つのインデックスをもつ（階数３の）テンソルと見なす．それを図 1(a) のように図示する．この表現ではノードと枝はそれぞれテンソ
ルとそのインデックスを表している．すると，MPS での行列積は２つの（階数３の）テンソルの１つのインデックスに対するテンソル縮約である：[M (s2 )M (s3 )]t1 t3 =!

M (s2 )t1 t2 M (s3 )t2 t3 .

(6)

t2

そして，行列積の結果として（階数４の）テンソルが得られる．これを図示したのが図 1(b) である．テンソル縮約をとるインデックス (t2 ) は各ノード（テンソル）間を結ぶ枝で表現されている．このような表現をつかって，MPS は図 1(c) のようにノードと枝からなるテンソルのネットワークとして，つまり，テンソルネッ
トワークとして表現される．

2(a)(b)[M (s)]ij ≡i M j[M (s2 )M (s3 )]t1 t3 =st1 M!

M (s2 )t1 t2 M (s3 )t2 t3t2

≡t2

s2(c)

s3t6T

≡

s1 s2 s3 s4 s5 s6

図1M t3M

t1

M

s1t2s2Mt3

M

s3t4s4

M

s5t5M

s6

MPS のテンソルネットワーク表現．(c) ６スピン系の MPS の確率振幅．

DMRG は，図 1(c) の１次元的なテンソルネットワークで表現される量子状態の中で最小エネルギーを探索する計算手法で，テンソルネットワークを用いた変分法であるといえる：min⟨Ψ|H|Ψ⟩,M (↑),M (↓)(7)

ここでは最適化するパラメータはテンソル M である．
テンソルネットワークを用いた変分法は量子状態のクラスをテンソルネットワークで定義するのが特徴である．従って，系の基底状態が仮定したテンソルネットワークで表現できるかどうかが非常に重要である．近年，量子情報の研究が進展し，テンソルネットワークのもつ表現能力に関する理解が進んだ．その中で，特に重要なのは量子的相関に関する表現能力である．古典的な系では系が一つの状態にある時，その状態は常に場所ごとの状態の直積になっている．つまり，理想的な状況ではある場所での測定が他の場所での状態に影響を及ぼさない．量子状態でも，例えば，次の重ね合わせの状態は直積状態といえる．

.
/1111√ | ↑↓ ⟩ − √ | ↑↑ ⟩ = | ↑ ⟩ ⊗ √ | ↓ ⟩ − √ | ↑ ⟩ .
2222(8)このような場合は量子的な相関は２スピン間にないといわれる．一方，例えば，次の状態は直積で表現できず，一つのスピンに対する測定は直ちに残りのスピンの状態に影響を及ぼす．

11√ | ↑↓ ⟩ − √ | ↓↑ ⟩.
22(9)このように古典系にはないような仕方で相関をもつ場合，２スピン間には量子的な相関があるといわれる．また，量子的な相関をエンタングルメントと呼ばれている．
量子系を２分割したとき，部分系間のエンタングルメントを定量的に定義するものがエンタングルメントエントロピーである．いま，系を A と B の２つの部分系に分割する．部分系 A の縮約密度演算子 ρA は，全系の密度演算子 ρ の部分系 B に関する部分トレースを取ったもので定義される．部分トレースはトレースを一般化したもので，二つの系 A と B の合成系の演算子に対し，B に対する部分トレースは次のように定義される．

TrB [|sA ⟩⟨s′A | ⊗ |sB ⟩⟨s′B |] ≡ (⟨sB |s′B ⟩)|sA ⟩⟨s′A |.
3(10)今の場合，基底状態 |Ψ⟩ の密度演算子は ρ = |Ψ⟩⟨Ψ| である．従って，基底状態が A と B の直積状態|Ψ⟩ ≡ |ΨA ⟩ ⊗ |ΨB ⟩ であれば，ρA = |ΨA ⟩⟨ΨA | となり，固有値１の固有状態が１つだけ存在するが，A と Bがエンタングルメントしていれば複数の固有状態を持つことがわかる．このような縮約密度行列の固有状態の分布を次のエントロピー的な量として定義したものがエンタングルメントエントロピー SA （定義より常にSA = SB ）である．
SA = −Tr[ρA log ρA ].

(11)

エンタングルメントエントロピーの最も重要な法則は，部分系 A と B の境界の大きさに比例するというエリア則である．一方，テンソルネットワークで表現される量子状態のエンタングルメントエントロピーの上限値は，部分系 A と B を切り離す最小カット枝数に比例するということが知られている．従って，エンタングルメントエントロピーのエリア則を満たすためには，枝が境界の広さに比例して必要になる．MPS が１次元系をうまく表現できたのは１次元では境界は多くても２箇所であり定数だからである．一方，２次元以上の量子系では，境界が系の断面の大きさを持つため，MPS をそのまま適用できない．実際，図 2(a) のように２次元に MPS を無理に適用すれば，エリア則を満たさないため，その表現能力は非常に制限される．この欠点を克服するために近接点格子間を新しい枝で結んだネットワークが PEPS と呼ばれるテンソルネットワークでMPS を２次元以上の系に拡張したものとなっている (図 2(b))．PEPS はエリア則を満たすため，２次元以上の量子状態をコンパクトに記述できる可能性をもつとして注目を集めている．

(a)(b)

A
L

図2χ

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

M

MMPS(a) と PEPS(b) の最小カット. (b) 部分系間の枝数が境界の長さに比例．

PEPS は MPS の拡張であるが，ツリーネットワークの拡張として MERA が提案されている．ツリーネッ
トワークは実空間繰り込み的な構造をもつテンソルネットワークである．例えば，図 3(a) は，３サイトを１サイトに繰り込むテンソル（下に３本，上に１本の枝がある三角形テンソル）を積み重ねて構成されている．
２次元以上でも同様であるが，ツリーネットワークでは，最小カット数を１にできる場合がある（図 3(a)）．
従って，エリア則を満たさない．一方，MERA は繰り込む前に周辺領域との間に Disentangler と呼ばれるテンソルを挿入し，領域間のエンタングルメントを局所的に弱めてから実空間繰り込みを行う：例えば，図 3(b)では四角形テンソルが Disentangler．常に繰り込む領域間に挿入するので，最小カット枝数は境界の大きさに比例することになる．従って，PEPS と同様に MERA もエリア則を満たす．
以上のテンソルネットワーク PEPS と MERA は共にエリア則を満たすために，高次元量子系の基底状態計算に対するブレークスルーをもたらすものとして非常に注目を集めている．

4

(a)(b)A

図3A

ツリーネットワーク (a) と MERA(b) の最小カット．赤い枝は部分系 A を切り離すための最小カット枝．

2.2 古典系の自由エネルギー
古典多体系の振る舞いを研究する上で最も重要な量は自由エネルギーである．自由エネルギーは全ての物理量の母関数である．しかし，自由エネルギーはボルツマン重みの全状態についての和であるので，状態数がシステムサイズに対して指数関数的に増大するために，大きな系の厳密な計算は難しい．

exp(−F/T ) =!

exp[−H(S)/T ].

(12)

S

近似的な数値計算手法としては，量子系と同じくモンテカルロ法が用いられることがあるが，最近，テンソルネットワークを用いて計算する手法が発展しつつある．
一般的に物理系では全エネルギーは局所的な相互作用エネルギーの和になっている．例えば，正方格子点上の小さな磁性体の集合系を磁性体の古典モデルとして考える．格子点 i 上の小磁性体の状態を si とすると，全エネルギーは例えば隣接格子点対 ⟨ij⟩ に対する局所相互作用エネルギー H(si , sj ) を用いて次のように定義される．

H=!

H(si , sj ).

(13)

⟨ij⟩

格子点を図 4(a) にあるように，４点づつのグループに分割し，各４点に対して次の階数４のテンソルを定義する：

0As1 s2 s3 s4 = exp −!

1H(si , s(i mod 4)+1 )/T .

i(14)すると，自由エネルギーはテンソル As1 s2 s3 s4 らのテンソル縮約で書き直すことができ，図 4(b) のように 45度傾いた PEPS になる．これ以外にも，古典形の自由エネルギーを PEPS に射影する方法は複数存在し，常に古典 d 次元系の自由エネルギーは d 次元 PEPS に変換できることが知られている [4]．

3 テンソルネットワークの計算手法
多体問題のテンソルネットワークによる再定式化をここまででみてきた．残された問題はテンソルネットワークによって記述された計算を行う手法である．本章では，MERA と PEPS という２つの代表的なテンソルネットワークによって記述される計算問題の計算手法を紹介する．

5(a)(b)AijAs1AAexp(−F/T ) =As4A

AAAs2s3A

図 4 古典系の自由エネルギーの PEPS 表現．

3.1 MERA でのテンソル縮約
MERA は臨界点など幅広い状態を記述できるため，量子系の基底状態計算への応用がいくつか始まっている．以下では主に三角格子上に定義された系を考える．
三角格子上の量子スピンモデルは BEDT-TTF系や Pd(dmit)2 系等の有機モット絶縁体の理論

χ3モデルとして研究されてきた．一般的に２次元以上の系に有効であると考えられている量子モンテ

wカルロ法を形式的に適用すると，サンプルの重みχ2が負になるいわゆる負符号問題が存在し，深刻なv

桁落ちによる精度の悪化が起きる．他にも三角格子のように幾何学的なフラストレーションをもつ

χ1

量子系の理論モデルはさまざまなものが存在する

uが，一般的に量子モンテカルロ法の負符号問題が存在するため，それらの数値的研究は停滞してきた．しかし，前章でみたようにテンソルネットワー
クを用いた変分法はそのような問題がないため，DMRG のように強力な手法になるのではないか

図5と期待されている．例えば，MERA を用いた三角三角格子のエンタングルメント繰り込み．

格子上の量子スピンモデルの研究 [5] では図 5 の３次元テンソルネットワークが用いられている．図 5 は 19 格子点を 1 格子点に繰り込む部分だけを抜き出してきたもので，全部で３層からなっており第１層と第２層は Disentangler 層となっている．このような局所的なエンタングルメントを弱めてから実空間繰り込みを行うことを，エンタングルメント繰り込みと呼ぶ．
テンソルネットワークを用いた変分法では，テンソルネットワークで記述される量子状態の物理量の期待値計算が主たる計算になる．例えば，基底状態計算のようにエネルギー最低状態を求める場合は，エネルギーの期待値計算が必要である．特に式 (4) から局所ハミルトニアンの期待値計算ができれば全エネルギー計算はできる．例えば，２格子点上の局所演算子等の期待値の MERA での計算を考える．

⟨A⟩ = ⟨ΨMERA |A|ΨMERA ⟩.

(15)

この計算は量子状態がテンソルネットワークで表されているので，式 (15) 自身もテンソルネットワークで表6現することができる．例えば，１次元系では図 6(b) のテンソルネットワークになる．更に，MERA では各テンソルは上下方向に Isometry と呼ばれる性質をもつものに限定するため，図 6(a) の簡単化が行え，図 6(b)の計算は図 6(c) のサブテンソルネットワークの計算に帰着される．このサブテンソルネットワークを causalcone (CC) と呼ぶ．CC の大きさはシステムサイズに比例せず，そのテンソル縮約の計算量はシステムサイズによらない．これは PEPS にはない MERA の特徴である．従って，CC の計算量は，テンソルのインデックスの自由度の多項式になる [6]．三角格子上の MERA の場合，例えば，図 7 のような CC の計算が必要にな
Top tensor(c)(a)(b)AAA

図6

MERA の期待値計算．

図7

三角格子 MERA の causal cone 計算．

る．CC のテンソル縮約計算の計算量は一般的に高次の多項式になり，この場合では 14 次にもなる．実際，三角格子の研究 [5] では，空間異方性のあるモデルを計算しており，192 × 6 = 2166 点からなるユニットセルをベースにした計算を行っている．ここから生成される１つの CC の計算では，テンソル縮約の途中に非常に大きなテンソルが中間生成されてしまうが，１ノードのメモリー上限（この例では 28 ギガ）を超えて計算を行うにはかなり複雑なコードを必要とする．ただ，各々の CC 計算は完全に独立なので並列化可能である．
そのため，現在の所，BLAS3 を用いた１ノード単体での CC コードを並列に駆動している．テンソルのインデックスの自由度の増大に対する CC 計算の計算（メモリー）量の急激な増大が現在の MERA 変分法の弱点であり，モンテカルロサンプリングによる統計的評価法も提案されている [7]．
このような厳しい制限があるが，MERA を用いることにより，空間異方性のあるパラメータ領域では，スピンの向きが滑らかに回転していくスパイラル状態が幅広く安定的に存在することが示された．特に，スパイラルのねじれ角の MERA の結果（図 8）は，他の解析的手法（級数展開）と良い一致を示し，スパイラル相の存在をサポートする強力な数値的証拠を示している．
これ以外にも，カゴメ格子 [8] やシャストリー・サーザーランド格子 [9] 等，従来法では計算が困難であった格子系でも MERA を用いた研究が行われており，今後の更なる応用が期待されている．

7

図 8 空間異方性のある三角格子反強磁性ハイゼンベルグモデルの磁性モーメント間のねじれ角．横軸は鎖内と鎖間の結合定数（J1 と J2 ）の比で空間異方性を表す．縦軸は J2 軸方向の磁性モーメント間のねじれ角．丸と三角が 2166 格子点からなるユニットセルを用いた MERA の結果で，実線が級数展開法の結果．

3.2 PEPS 計算での HOSVD の利用
古典系の自由エネルギーは 2.2 節で述べたように PEPS を用いて表現される．従って，PEPS 計算を行うことで自由エネルギーを計算することはできるが，PEPS の厳密な計算のための計算量は，d 次元 PEPS の場
合，システムサイズ L に対して指数 Ld−1 程度の指数関数的な増大を示す．従って，様々な近似的なテンソル縮約法が提案されている．本節では特に最近提案された HOSVD をもちいた計算手法について紹介する [10]．

HOSVD は行列の特異値分解のテンソルへの拡張になっている．テンソル T を次のような形のテンソル積に分解する (Tucker 分解)．

Tijkl =!

(1)

(2)

(3)

(4)

Si′ j ′ k′ l′ Uii′ Ujj ′ Ukk′ Ull′ ,(16)

i′ j ′ k ′ l ′ここでは，階数４のテンソルを例として分解している．U (t) はユニタリー行列で，S はコアテンソルと呼ばれている．コアテンソルについては，次の関係式が成立している．

⟨Sin =p , Sin =q ⟩ = δpq||Sin =1 || ≥ ||Sin =2 || ≥ · · · ,||A|| ≡2(17)⟨A, A⟩,(18)

ここで，Sik =p は，k 番目のインデックスの値が p である部分コアテンソルであり，部分コアテンソル同士の内積は，⟨A, B⟩ =3ij··· Aij··· × Bij··· で定義される．階数が２の場合は通常の特異値分解になる．それ以外′の場合でも，式 (18) より，中間変数の i などを打ち切れば，元のテンソル T の良い近似表現を構成できる．
この Tucker 分解をテンソルネットワークとして図示したものが図 9(a) である．
正方格子 PEPS の HOSVD を用いた近似的計算では，まず最初に横方向に隣接する２テンソルのテンソル積を１テンソルに置き換える．さらに縦方向にも同様な置き換えを行う．これでも元の正方格子の縦横半分の正方格子 PEPS が出来上がる．この手続きを正方格子のサイズが十分小さくなるまで繰り返す．この時，最も重要なポイントは２テンソルのテンソル積の結果を１テンソルに近似する所である．HOSVD はこの部分

8

(a)y′y′
V'xx′

TxU

SU'x′

Vyy(b)
y1′y2′

M

My′y1′ y2′y′V'xx′xM'x′xU

SU'x′xM''x′

Vy1y2y1 y2
図9yyHOSVD による PEPS 計算．

で用いられる．具体的には，図 9(b) にあるように，HOSVD で求めたユニタリー行列のインデックスを適当な所で打ち切ることで元テンソル積の近似を行う．

HOSVD を用いた手法のメリットは，従来のテンソル繰り込み群 (Tensor renormalization group, TRG)
よりも高精度である点と，任意次元の PEPS に適用できる点である．代表的な古典多体系である古典 Ising モ
デルでの TRG との比較は文献 [10] に詳しく紹介されている．３次元 PEPS に適用可能であるということは，有限温度の２次元量子系などの計算も直接おこなることを意味し，今後の展開が期待されている．

4 まとめと今後の展望
本稿では，物性研究で取り扱われる多体問題に対して，テンソルネットワークを用いた計算手法がどのように適用されるかをいくつかの例を通して紹介した．従来法では取り扱いが困難であった重要な問題がテンソルネットワークをベースにした計算手法を武器に現在活発に研究されている．これらの手法を用いた新しい物性研究は，計算手法自身や応用研究共に非常に活発に進められており，今後のブレークスルーが期待されている新しい研究分野を作りつつある．
ただ，BLAS や LAPACK や ScaLAPACK 等の高性能の線形計算ライブラリーの整備に比して，テンソル計算の高性能ライブラリーはまだ十分とはいえない．特に大規模な並列計算機向けの高性能ライブラリーの整備が望まれている．更に，HOSVD などテンソル分解に関する研究もテンソルネットワークの今後の進展に大きな影響を与えると考えられる．計算物理と応用数学のコラボレーションがこのような現状を打破し，従来研究を超えた新しい成果が得られることが望まれている．

参考文献[1] G Vidal. Entanglement Renormalization. Physical Review Letters, 99(22):220405, November 2007.
[2] F. Verstraete, M Wolf, D Perez-Garcia, and J Cirac. Criticality, the Area Law, and the Computational Power of Projected Entangled Pair States. Physical Review Letters, 96(22):220601, June2006.
[3] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. A Multilinear Singular Value Decom-9position. SIAM Journal on Matrix Analysis and Applications, 21(4):1253–1278, January 2000.
[4] H. H. Zhao, Z Y Xie, Q. N. Chen, Z. C. Wei, J W Cai, and T. Xiang. Renormalization of tensornetwork states. Physical Review B, 81(17), May 2010.
[5] Kenji Harada. Numerical study of incommensurability of the spiral state on spin-1/2 spatiallyanisotropic triangular antiferromagnets using entanglement renormalization. Physical Review B,
86(18):184421, November 2012.
[6] G Evenbly and G Vidal.

Algorithms for entanglement renormalization.

Physical Review B,

79(14):144108, April 2009.
[7] Andrew J Ferris and Guifre Vidal. Variational Monte Carlo with the multiscale entanglement
renormalization ansatz. Physical Review B, 85(16):165146, April 2012.
[8] G Evenbly and G Vidal.

Frustrated Antiferromagnets with Entanglement Renormalization:

Ground State of the Spin-1/2 Heisenberg Model on a Kagome Lattice. Physical Review Letters,
104(18):187203, May 2010.
[9] J Lou, T Suzuki, K Harada, and N Kawashima. Study of the Shastry Sutherland Model Using
Multi-scale Entanglement Renormalization Ansatz. arXiv:1212.1999.
[10] Z Xie, J. Chen, M Qin, J Zhu, L Yang, and T. Xiang. Coarse-graining renormalization by higherorder singular value decomposition. Physical Review B, 86(4):045139, July 2012.

10