6. アルゴリズム

6.1. テンソルネットワーク状態

テンソルネットワーク状態 (Tensor network states (TNS)) とは小さなテンソルの積、繋がりで表現された変分波動関数です [TNS] 。例えば、\(N\) 個の \(S=1/2\) 量子スピン系では、その波動関数は直積状態の基底を用いて、

\[|\Psi\rangle = \sum_{s_i = \uparrow,\downarrow} \Psi_{s_1,s_2,\dots,s_N} |s_1,s_2,\dots,s_N\rangle\]

と表せます。テンソルネットワーク状態では、この展開係数 \(\Psi_{s_1,s_2,\dots,s_N}\) はテンソルのネットワークで表現され、例えば

\[\Psi_{s_1,s_2,\dots,s_N} = \mathrm{tTr}\left[T^{(1)}[s_1]T^{(2)}[s_2]\cdots T^{(N)}[s_N]\right],\]

と書くことができます。ここで、 \(\mathrm{tTr}[\dots]\) はテンソルネットワークの縮約を表し、 \(T^{(i)}[s_i]\) はテンソルを表しています。 行列積状態 (matrix product state (MPS)) と呼ばれるテンソルネットワーク状態の場合 [MPS] , \(T^{(i)}[s_i]\)\(s_i\) が与えられると行列になっていて、 \(\mathrm{tTr}[\dots]\) はこの場合、通常の行列積になります:

\[\Psi_{s_1,s_2,\dots,s_N}^{\mathrm{MPS}} = T^{(1)}[s_1]T^{(2)}[s_2]\cdots T^{(N)}[s_N],\]

ここで、 \(T^{(1)}[s_1]\)\(T^{(i)}[s_i] (i\neq 1, N)\)\(T^{(N)}[s_N]\) の形状はそれぞれ、 \(1 \times D_1\)\(D_{i-1} \times D_{i}\)\(D_{N-1} \times 1\) であると仮定しました。

このようなTNSを基底状態波動関数の近似に用いる場合、その精度は \(D_i\) によって制限されます。 \(D_i\)ボンド次元 と呼ばれています。テンソルネットワークのダイアグラム表記を用いると、 MPSは

../_images/MPS.png

のように描くことができます。このMPSは有限系の波動関数を表していますが、同様にして、無限系の波動関数を表わす、無限に長いMPSを考えることもできます。特に、波動関数が格子の(ある長さの)並進に対して対称性を持っている場合、少数の独立なテンソルを繰り返すことで、無限系のMPS (infinite MPS (iMPS)) を作ることができます。 2サイトの並進対称性の場合、このような iMPS は

../_images/iMPS.png

と表せます。ここで、同じ色のテンソルは、同じ要素を持つテンソルです。

TeNeS では、2次元の無限に広がったテンソル積状態 (infinite tensor product states (iTNS)) を取り扱います。この iTPS は iMPS の高次元への自然な拡張になっています。 TeNeSでは、並進対称性をもった正方格子テンソルネットワークを仮定しており、ダイヤグラムでは、

../_images/iTPS.png

のように描けます。TeNeSではこのiTPSを用いて、2次元量子多体系の基底状態を近似的に計算します。なお、正方格子テンソルネットワークは、正方格子模型だけでなく、適切なマッピングにより、ハニカム格子模型、三角格子模型など様々な2次元格子模型に適用できます。

6.2. iTPS の縮約

あるTNSが与えられた時に、そのTNSでの期待値、 \(\langle \Psi|O|\Psi\rangle/\langle \Psi|\Psi\rangle\) を計算するためには、 一般に、 \(\langle \Psi|O|\Psi\rangle\)\(\langle \Psi|\Psi\rangle\) という二つの量に対応するテンソルネットワークの縮約計算が必要になります。 例えば、 \(\langle \Psi|\Psi\rangle\) に対応するテンソルネットワークは

../_images/iTPS_braket.png

で与えられます。この形のテンソルネットワークは、しばしば、ダブルレイヤー(double layered)テンソルネットワークと呼ばれます。 ダブルレイヤーテンソルネットワークの縮約計算は、通常、非常に大きな計算コストが必要です。MPS(や iMPS)の場合には、幸いにも、局所的なテンソルで構成される転送行列を考えることなどによって、効率的に計算することができます。しかし、TPS(や iTPS)の場合、厳密な縮約計算は小さいクラスター(又は小さい半径の無限シリンダー)を除いてほぼ不可能で、通常、近似的な縮約計算法を用います。TeNeSでは、角転送行列繰り込み群法(corner transfer matrix renormalization group (CTMRG) [CTMRG] と呼ばれる、無限に広がったダブルレイヤーテンソルネットワークを 角転送行列エッジテンソル を用いて近似する方法を採用しています。

ダブルレイヤーテンソルネットワークを局所的に縮約したテンソル

../_images/double_tensor.png

を使って単純化すると、角転送行列表現に対応するテンソルネットワークダイアグラムは、

../_images/CTM.png

と表されます。角転送行列とエッジテンソルは、

../_images/CandE.png

のように定義されています。角転送行列表現の精度は、ダイアグラム中で太線で表現した、角転送行列のボンド次元 \(\chi\) によって制限されます。

CTMRGのアルゴリズムでは、 角転送行列とエッジテンソルに局所的なテンソルを 吸収 していくことでそれらをアップデートし、結果が収束するまで繰り返します。例えば、 left move と呼ばれる吸収手続きは、ダイアグラムでは

../_images/LeftMove.png

と表されます。このダイアグラムに現れる プロジェクター は、いくつかの方法で計算することができ [CTMRG] 、自由度を \(\chi\) に減らす働きをします。

ボンド次元 \(D\) のiTPSを用いて、ボンド次元 \(\chi\) の角転送行列表現を考える場合、CTMRGの計算コストは、 \(O(\chi^2 D^6)\)\(O(\chi^3 D^4)\) の大きな方でスケールします。 ここで、ダブルレイヤーテンソルネットワークのボンド次元は、局所縮約したテンソルを用いる表現では、 \(D^2\) になっていることに注意してください。このため、通常、 \(\chi\)\(\chi \propto O(D^2)\) のように \(D^2\) に比例して増やします。この条件では、CTMRGの計算コストは \(O(D^{10})\) になり、メモリ量は \(O(D^{8})\) になります。 なお、ここで述べた計算コストを得るためには、疎行列の特異値分解(SVD)を用いる必要があります、代わりに、密行列のSVDを用いる場合、計算コストは \(O(D^{12})\) となります。

いったん収束した角転送行列とエッジテンソルを得れば、 \(\langle \Psi|O|\Psi\rangle\) も効率的に計算することができます。例えば、局所磁化 \(\langle \Psi|S^z_i|\Psi\rangle\) は、

../_images/Sz.png

のように表わされ、同様に最近接相関 \(\langle \Psi|S^z_iS^z_{i+1}|\Psi\rangle\)

../_images/SzSz.png

と表現することができます。また、ダイアグラムの2番目の表記を用いることで、任意の2サイト演算子の期待値も計算できることがわかります。このようなダイグラムを任意の演算子に対して描くことは可能ですが、クラスターが大きくなるとその縮約計算に必要となる計算コストが莫大になることに注意してください。

6.3. iTPSの最適化

iTPSを基底状態の変分波動関数として用いる場合、iTPSが最小のエネルギー期待値

\[E = \frac{\langle \Psi|\mathcal{H}|\Psi\rangle}{\langle \Psi|\Psi\rangle},\]

を与えるように、テンソルを最適化する必要があります。ここで、 \(\mathcal{H}\) は対象系のハミルトニアンを表しています。 TeNeSでは、虚時間発展(the imaginary evolution (ITE))法と変分最適化(the variational optimization)法という二つの手法のうち、前者の ITE を採用しています。TeNeS では、iTPSの範囲での近似的な虚時間発展

\[|\Psi^{\mathrm{iTPS}} \rangle \simeq e^{-T \mathcal{H}} |\Psi_0\rangle,\]

を考えます。ここで、 \(|\Psi_0 \rangle\) は任意の初期 iTPS です。. もし、 \(T\) が十分に大きければ、左辺の \(|\Psi^{\mathrm{iTPS}}\rangle\) は基底状態の良い近似になっていると考えることができます。

TeNeSでは、ハミルトニアンは短距離の二体相互作用の和で

\[\mathcal{H} = \sum_{\{(i,j)\}}H_{ij},\]

のように表されていると仮定し、小さな時間刻み \(\tau\) の虚時間発展演算子に対してSuzuki-Trotter 分解

\[e^{-\tau \mathcal{H}} = \prod_{\{(i,j)\}} e^{-\tau H_{ij}} + O(\tau^2).\]

を適用します。ここでは、一次の近似を考えましたが、より高次の分解を考えることもできます。Suzuki-Trotter 分解の形を用いることで、虚時間発展は

\[e^{-T \mathcal{H}} |\Psi_0\rangle = \left( \prod_{\{(i,j)\}} e^{-\tau H_{ij}}\right)^{N_{\tau}} |\Psi_0\rangle + O(\tau),\]

のように書き下すことができます。ここで、 \(N_{\tau} = T/\tau\) は十分に小さな \(\tau\) での虚時間発展のステップ数です。 この式の右辺を計算するために、 \(\prod_{\{(i,j)\}}\) の積をいくつかの部分集合に分解します。それぞれの部分集合内では, (局所的な)虚時間発展演算子はお互いに交換し、考えているiTPSと同じ並進対称性を持っているとします。例えば、2サイトの iMPS で、1元系の最近接相互作用ハミルトニアンを考えた場合、二つの部分集合を用いて、

../_images/iMPS_ITE.png

のように虚時間発展を分解することができます。

次に、それぞれの虚時間発展演算子の部分集合を適用した波動関数を、ボンド次元 \(D\): の新しいiTPSとして

\[|\Psi_{\tau}^{\mathrm{iTPS}} \rangle \simeq \prod_{\{(i,j) \in \mathrm{subset}_n \}}e^{-\tau H_{ij}} |\Psi^{\mathrm{iTPS}}\rangle,\]

のように近似します。ここで \(\prod_{\{(i,j) \in \mathrm{subset}_n \}}\)\(n\) 番目の部分集合ないの演算子の積を表し、 \(|\Psi_{\tau}^{\mathrm{iTPS}}\rangle\) は新しい iTPS です。ダイアグラムを用いるとこの式は、

../_images/iMPS_ITE_iMPS.png

のように表現できます。一般に、 \(e^{-\tau H_{ij}}\) をかけることで 厳密な iTPS 表現のボンド次元は増大してしまうことに注意してください。したがって、虚時間発展のシミュレーションを安定して継続するためには、ボンド次元をある一定値 \(D\) まで毎回打ち切る ( truncate ) 必要があります。

素朴には、効率的な打ち切りは、最小化問題

\[\min \left \Vert |\Psi_{\tau}^{\mathrm{iTPS}} \rangle -\prod_{\{(i,j) \in \mathrm{subset}_n \}} e^{-\tau H_{ij}} |\Psi^{\mathrm{iTPS}}\rangle \right \Vert^2.\]

を解くことで行えます。しかし、この最小化問題を解く計算コストは、主にiTPSの並進対称性で問題が非線形問題になっているために、非常に膨大になってしまいます。そこで、通常は、代わりの問題として、局所的な一つの虚時間発展演算子だけを適用して、それを近似する iTPS \(|\Psi_{\tau}^{\mathrm{iTPS}}\rangle\) を探す問題を考えます。 ここで、新しいiTPSでは、元の \(|\Psi^{\mathrm{iTPS}}\rangle\) と比較して、数個のテンソルだけが変更されています。 この局所的な最小化問題は

\[\min \left \Vert |\Psi_{\tau}^{\mathrm{iTPS}} \rangle - e^{-\tau H_{ij}} |\Psi^{\mathrm{iTPS}}\rangle \right \Vert^2\]

と書くことができます。一次元の最近接相互作用の場合、この最小化問題に対応するダイアグラムは、

../_images/iMPS_ITE_local.png

で与えられます。

差の二乗ノルム \(\left \Vert |\Psi_{\tau}^{\mathrm{iTPS}} \rangle - e^{-\tau H_{ij}} |\Psi^{\mathrm{iTPS}}\rangle \right \Vert^2\) は、例えば CTMRG 等を使うことで効率的に計算できるため、この最適化問題は簡単に解くことができます [ITE] 。ここで新しく得られる iTPS は並進対称性を破っていますが、アップデートされたテンソルを他の場所に コピー することで、並進対称な iTPS を作ることができます。

../_images/Copy.png

このiTPSは元の最小化問題の近似解だと考えることができます。このような虚時間発展の方法は、 full update 法と呼ばれます。 full update 法の計算の大部分は CTMRG であり、SVDの方法に応じて、計算コストは \(O(D^{10})\) または \(O(D^{12})\) でスケールします。

Simple update 法は虚時間発展を用いた、より計算コストの小さい最適化手法です。Simple update法では、CTMRGによる重い計算を避けるために, 波動関数全体ではなく、局所的なテンソルネットワークを考えます [SimpleUpdate] 。例えば 最近接相互作用の場合には、以下のような局所的な最適化問題を考えます。

../_images/Simple_opt.png

このダイアグラムでは、 \(\lambda_i\) は非負の対角行列を表していて、これはボンド \(i\) の先にある無視した環境を表わす平均場だと考えることができます。 \(\lambda_i\) の具体的な定義は後で与えられます。 このダイアグラムが表わす最適化問題は、テンソル二つと虚時間発展演算子一つが一体となった行列の低ランク近似と見做すことができるため、SVDを用いて解くことができます。この手続きは、ダイアグラムを用いて、

../_images/Simple_update.png

と表すことができます。計算途中のSVDで出てきた行列の特異値は、次のステップでの平均場 \(\lambda\) として利用されます。Simple update法の計算コストは、行列を構成する前にQR分解を行うことで、 \(O(D^{5})\) になります [QR] 。したがって、simple update法はfull update法よりもずっと計算コストが軽くなっています。

ただし、simple update法はfull updateよりも計算コストが小さいですが、simple update法は初期状態依存性が強く、また、最終結果の局所磁化の大きさを過剰評価する問題が知られています。したがって、未知の問題に適用する場合には、得られた結果を慎重に検証する必要があります。

参考文献

[TNS] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals. of Physics 349, 117 (2014). link; R. Orús, Tensor networks for complex quantum systems, Nature Review Physics 1, 538 (2019). link; 西野友年、大久保毅 テンソルネットワーク形式の進展と応用, 日本物理学会誌 72, 702 (2017). link; 大久保毅 テンソルネットワークによる情報圧縮とフラストレート磁性体への応用, 物性研究・電子版 7, 072209 (2018) link.

[MPS] U. Schollwcök, The density-matrix renormalization group in the age of matrix product states, Annals. of Physics 326, 96 (2011). link

[CTMRG] T. Nishino and K. Okunishi, Corner Transfer Matrix Renormalization Group Method, J. Phys. Soc. Jpn. 65, 891 (1996).; R. Orús and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009). link ; P. Corboz et al., Competing States in the t-J Model: Uniform d-Wave State versus Stripe State, Phys. Rev. Lett. 113, 046402 (2014). link

[ITE] J. Jordan et al., Classical Simulation of Infinite-Size Quantum Lattice Systems in Two Spatial Dimensions, Phys. Rev. Lett. 101, 250602, (2008). link; R. Orús and G. Vidal, Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction, Phys. Rev. B 80, 094403 (2009). link

[SimpleUpdate] H. G. Jiang et al., Accurate Determination of Tensor Network State of Quantum Lattice Models in Two Dimensions, Phys. Rev. Lett. 101, 090603 (2008). link

[QR] L. Wang et al., Monte Carlo simulation with tensor network states, Phys. Rev. B 83, 134421 (2011). link