🙆‍♀️

テンソルネットワーク入門

2023/09/17に公開

本記事は、初学者向けにテンソルネットワークについて導入したものである。
コードは、ITensor.jlのソフトウェア論文The ITensor software library for tensor network calculationsを参照した。

テンソルとは？？

テンソルとは、ある連続関数を離散化したもので、それを多次元配列で表現したものである。
例えば、1階のテンソルならば、1変数の連続関数の空間を離散化し、グリッド上で評価した関数である。これはベクトルである
2階のテンソルならば、2変数関数の空間を離散化し、そのグリッド点上での値である。これは行列である。
同様に3階,,と一般化できる。
高階のテンソルになればなるほど、空間の次元が上がるので、それに伴いグリッド点の数は、指数関数的に増大することに注意する。

ITensor.jlを使っていく前に、もし機能で困ったら、以下のhelp機能を使うということを覚えておく（重要）

julia> using ITensors
julia> ?
help?> Index

テンソルをダイアグラムでかくと、

さて、ITenosr.jlを使って、テンソルを表現していこう。

using ITensors
i = Index(4)
j = Index(4)
k = Index(4)

S = ITensor(1.) # scalar
A = ITensor(i) 
B = ITensor(i,j)
C = ITensor(i,j,k)

試しに行列をITensorで作ってみる

i = Index(2,"index_i")
j = Index(2,"index_j")

M = [1. 2;
     3 4]
T = ITensor(M, i, j)

テンソルの要素指定は、

T[i => 1, j => 1] = 3.3
M[1, 1] == 3.3
T[i => 1, j => 1] == 3.3

#別の例
C = ITensor(i,j,k)
C[i=>2,j=>3,k=>1] = 0.837

ランダムテンソルは、

#random tenosor
T = randomITensor(i,j,k)

onehot

i = Index(2,"i")
A = onehot(i=>2)
# A[i=>2] == 1, all other elements zero

# Specify the element type
A = onehot(Float32, i=>2)

j = Index(3,"j")
B = onehot(i=>1,j=>3)
# B[i=>1,j=>3] == 1, all other element zero

テンソルのみかた

高階のテンソルは、複数のlocal index(足）をまとめることで、ある行列としてみることができる。
この考え方は、のちのMPSを作るところで出てくる。

C = combiner(i,k; tags="c")
T = randomITensor(i,j,k)
CT = C * T
@show CT

図は論文から引用。https://scipost.org/SciPostPhysCodeb.4/pdf

縮約

テンソルに対する演算に、縮約というものがある。これは、異なるテンソル間で共通のindex（足）に対して、和をとる操作である。

x = ITensor(j)
A = ITensor(i,j)
D = x * A

行列の掛け算は、

A = ITensor(i,j)
B = ITensor(j,k)

C = A * B

これは以下と等しい

C = A * transpose(B)

テンソルの和

A = randomITensor(i,j,k)
B = randomITensor(k,i,j)
C = A + B

定数倍をかけるには、

D = 4*A - B/2
F = A + 3.0im * B

テンソル積

i = Index(3)
j = Index(3)
k = Index(3)
α = Index(3)
β = Index(3)

A = ITensor(i, α)
B = ITensor(α, j, β)
C = ITensor(β, k)
T = A * B * C

#delta関数とか、面白い例はたくさんある。

テンソル分解

特異値分解を使ったテンソル分解は、

i = Index(3)
j = Index(3)
A = randomITensor(i,j)
U,S,V = svd(A,(i))
U*S*V ≈ A # true

cutoffは閾値で、maxdimは最大のボンド次元である。

U,S,V = svd(W,(i);cutoff=1E-8,maxdim=10)

切り捨て誤差は、切り捨てられた特異値の二乗の合計を、全特異値の二乗の合計で割った値として計算される。この値がεより小さくなるように特異値を切り捨てられる設定をしている。
それと同時に、保持される特異値の最大数が10以下であることとしている。

Matrix Product State (MPS)

ここでは、行列積状態（Matrix Product States）の概要を簡単に説明する。
$A$ を $n_1 \times n_2 \times \cdots \times n_L$ のテンソルとし、 $L$ 個のインデックス $i_1, i_2,\cdots, i_L$ を持つオブジェクトとする。これを $A(i_1, i_2, \cdots, i_L)$ と表記する。
MPS（Matrix Product State）は、単一の $L$ 階テンソル $A$ を $L$ 個の3階のテンソル $A^{(1)},\cdots, A^{(L)}$ の縮約によって近似するものである。

A(i_1, i_2, \ldots, i_L) \simeq \sum_{\alpha_0=1}^{D_0} \cdots \sum_{\alpha_L=1}^{D_L} A_{i_1,\alpha_0\alpha_1}^{(1)} A_{i_2,\alpha_1\alpha_2}^{(2)} \cdots A_{i_L,\alpha_{L-1}\alpha_L}^{(L)} \equiv A_{i_1}^{(1)} \cdots A_{i_L}^{(L)}

ここで、 $A^{(l)}$ は補助的な $n_l × D_{l-1} × D_l$ のテンソルである。localまたは物理的なインデックス $i_1, \cdots, i_L$ に加えて、ダミーまたは仮想的なインデックス $\alpha_0, \cdots, \alpha_L$ を導入する。インデックス $\alpha_l$ はテンソル $A^{(l-1)}$ と $A^{(l)}$ の間の「ボンド」を形成し、 $D_l$ はボンド次元と呼ばれる（定義により、 $D_0=D_L=1$ である）。
MPSのボンド次元 $D$ は、構成要素の中で最大のボンド次元として定義され、 $D = \mathrm{max}_lD_l$ である。
物理的なインデックスが固定された状態で、ボンド方向にテンソル $A^{(l)}$ の縮約をとらせて行列積を形成することで、MPSの縮約表記 $A_{i_1}^{(1)} \cdots A_{i_L}^{(L)}$ を可能にする。

MPSは、基本的に、以下の図のように、高階のテンソルに対して、SVDを何度も繰り返すことで作ることができる。

ソフトウェア上で、MPSを作るために、何度もSVDを繰り返す操作をユーザー側はする必要がない（教育的には手を動かすことも大事だと思うので、実際に手を動かしてMPSを作ってみたい方は、こちらの画像の特異値分解の記事を参照あれ。）

cutoff = 1E-8
maxdim = 10
T = randomITensor(i,j,k,l,m)
M = MPS(T,(i,j,k,l,m);cutoff=cutoff,maxdim=maxdim)