2. iTEBD

By : Ke Hsu, Kai-Hsin Wu

Time evolution block decimation is one of the most simple and sucessful Tensor network method [Vid07]. The core concept of this algorithm is to use the imaginary time evolution to find the best variational ansatz, usually in terms of Matrix product state (MPS).

Here, we use a 1D transverse field Ising model (TFIM) as a simple example to show how to implement iTEBD algorithm in Cytnx and get the infinite system size (variational) ground state.

Consider the Hamiltonain of TFIM:

\[H = J\sum_{ij} \sigma^{z}_i\sigma^{z}_j - H_x\sum_i \sigma^{x}_i\]

where \(\sigma^{x,z}\) are the pauli matrices. The infinite size ground state can be represent by MPS as variational ansatz, where the virtual bonds dimension \(\chi\) effectively controls the number of variational parameters (shown as orange bonds), and the physical bonds dimension \(d\) is the real physical dimension (shown as the blue bonds. Here, for Ising spins \(d=2\)).

Because the system has translational invariant, thus it is legit to choose unit-cell consist with two sites, and the infinite system ground state can be represented with only two sites MPS with local tensors \(\Gamma_A\) and \(\Gamma_B\) associate with the schmit basis and \(\lambda_A\), \(\lambda_B\) are the diagonal matrices of Schmidt coefficients as shown in the following:

../_images/itebd_MPS.png

Let’s first create this two-site MPS wave function (here, we set virtual bond dimension \(\chi = 10\) as example)

  • In Python

 1import cytnx
 2
 3chi = 10
 4A = cytnx.UniTensor([cytnx.Bond(chi),cytnx.Bond(2),cytnx.Bond(chi)],labels=["a","0","b"]);
 5B = cytnx.UniTensor(A.bonds(),rowrank=1,labels=["c","1","d"]);
 6cytnx.random.normal_(B.get_block_(),0,0.2);
 7cytnx.random.normal_(A.get_block_(),0,0.2);
 8A.print_diagram()
 9B.print_diagram()
10#print(A)
11#print(B)
12
13la = cytnx.UniTensor([cytnx.Bond(chi),cytnx.Bond(chi)],labels=["b","c"],is_diag=True)
14lb = cytnx.UniTensor([cytnx.Bond(chi),cytnx.Bond(chi)],labels=["d","e"],is_diag=True)
15la.put_block(cytnx.ones(chi));
16lb.put_block(cytnx.ones(chi));
17la.print_diagram()
18lb.print_diagram()

Output >>

-----------------------
tensor Name :
tensor Rank : 3
block_form  : false
is_diag     : False
on device   : cytnx device: CPU
            -------------
           /             \
    -1 ____| 10        2 |____ 0
           |             |
           |          10 |____ -2
           \             /
            -------------
-----------------------
tensor Name :
tensor Rank : 3
block_form  : false
is_diag     : False
on device   : cytnx device: CPU
            -------------
           /             \
    -3 ____| 10        2 |____ 1
           |             |
           |          10 |____ -4
           \             /
            -------------
-----------------------
tensor Name :
tensor Rank : 2
block_form  : false
is_diag     : True
on device   : cytnx device: CPU
            -------------
           /             \
    -2 ____| 10       10 |____ -3
           \             /
            -------------
-----------------------
tensor Name :
tensor Rank : 2
block_form  : false
is_diag     : True
on device   : cytnx device: CPU
            -------------
           /             \
    -4 ____| 10       10 |____ -5
           \             /
            -------------

Here, we use random::normal_ to initialize the elements of UniTensor A and B with normal distribution as initial MPS wavefuncion. The la, lb are the weight matrix (schmit coefficients), hence only diagonal elements contains non-zero values. Thus, we set is_diag=True to only store diagonal entries. We then initialize the elements to be all one for this weight matrices.

Note

In general, there are other ways you can set-up a trial initial MPS wavefunction, as long as not all the elements are zero.

2.1. Imaginary time evolution

To optimize the MPS for the ground state wave function, in TEBD, we perform imaginary time evolution with Hamiltonian \(H\) with evolution operator \(e^{\tau H}\). The manybody Hamiltonian is then decomposed into local two-sites evolution operator (or sometimes also called gate in quantum computation language) via Trotter-Suzuki decomposition, where \(U = e^{\tau H} \approx e^{\delta \tau H_{a}}e^{\delta \tau H_{b}} \cdots = U_a U_b\), \(U_{a,b} = e^{\delta \tau H_{a,b}}\) are the local evolution operators with \(H_a\) and \(H_b\) are the local two sites operator:

\[H_{a,b} = J\sigma^{z}_{A,B}\sigma^{z}_{B,A} - \frac{H_x}{2}(\sigma^{x}_A + \sigma^{x}_B)\]

This is equivalent as acting theses two-site gates consecutively on the MPS, which in terms of tensor notation looks like following Figure(a):

../_images/itebd_upd.png

Since we represent this infinite system MPS using the translational invariant, the Figure(a) can be further simplified into two step. First, acting \(U_a\) as shown in Figure(1) then acting \(U_b\) as shown in Figure(2). This two procedures then repeat until the energy is converged.

Here, let’s construct this imaginary time evolution operator with parameter \(J=-1\), \(H_x = -0.3\) and (imaginary) time step \(\delta \tau = 0.1\)

  • In Python

 1J  = -1.0
 2Hx = -0.3
 3CvgCrit = 1.0e-10
 4dt = 0.1
 5
 6## Create onsite-Op
 7Sz = cytnx.physics.pauli("z").real()
 8Sx = cytnx.physics.pauli("x").real()
 9I = cytnx.eye(2)
10
11## Build Evolution Operator
12TFterm = cytnx.linalg.Kron(Sx,I) + cytnx.linalg.Kron(I,Sx)
13ZZterm = cytnx.linalg.Kron(Sz,Sz)
14
15H = Hx*TFterm + J*ZZterm
16del TFterm, ZZterm
17print(H)
18
19eH = cytnx.linalg.ExpH(H,-dt) ## or equivantly ExpH(-dt*H)
20eH.reshape_(2,2,2,2)
21print(eH)
22H.reshape_(2,2,2,2)
23
24eH = cytnx.UniTensor(eH)
25eH.print_diagram()
26print(eH)
27
28H = cytnx.UniTensor(H)
29H.print_diagram()

Output>>

Total elem: 4
type  : Double (Float64)
cytnx device: CPU
Shape : (2,2)
[[1.00000e+00 0.00000e+00 ]
 [0.00000e+00 -1.00000e+00 ]]


Total elem: 4
type  : Double (Float64)
cytnx device: CPU
Shape : (2,2)
[[0.00000e+00 1.00000e+00 ]
 [1.00000e+00 0.00000e+00 ]]


Total elem: 16
type  : Double (Float64)
cytnx device: CPU
Shape : (4,4)
[[-1.00000e+00 3.00000e-01 3.00000e-01 0.00000e+00 ]
 [3.00000e-01 1.00000e+00 0.00000e+00 3.00000e-01 ]
 [3.00000e-01 0.00000e+00 1.00000e+00 3.00000e-01 ]
 [0.00000e+00 3.00000e-01 3.00000e-01 -1.00000e+00 ]]

-----------------------
tensor Name :
tensor Rank : 4
block_form  : false
is_diag     : False
on device   : cytnx device: CPU
            -------------
           /             \
     0 ____| 2         2 |____ 2
           |             |
     1 ____| 2         2 |____ 3
           \             /
            -------------

Note

  1. Since \(U_a\) and \(U_b\) have the same content(matrix elements) but acting on different sites, we only need to define a single UniTensor.

  2. Here as a simple example, we directly convert a cytnx.Tensor to cytnx.UniTensor, which we don’t impose any bra-ket constrain (direction of bonds). In general, it is also possible to give bond direction (which we refering to tagged) that constrain the bonds to be more physical. See Github example/iTEBD/iTEBD_tag.py for demonstration.

  3. In general, the accurate ground state can be acquired with a higher order Trotter-Suzuki expansion, and with decreasing \(\delta \tau\) along the iteraction. (See [Vid07] for further details), Here, for demonstration, we use fixed value of \(\delta \tau\).

Tip

Here, physics.pauli returns complex type cytnx.Tensor. Since we know pauli-z and pauli-x should be real, we use .real() to get the real part.

2.2. Update procedure

Now we have prepared the initial trial wavefunction in terms of MPS with two sites unit cell and the time evolution operator, we are ready to use the aformentioned scheme to find the (variational) ground state MPS. At the beginning of each iteration, we evaluate the energy expectation value \(\langle \psi | H | \psi \rangle / \langle \psi | \psi \rangle\), and check the convergence, the network is straightforward:

../_images/itebd_contract.png ../_images/itebd_energy.png

  • In Python

 1A.relabels_(["a","0","b"])
 2B.relabels_(["c","1","d"])
 3la.relabels_(["b","c"])
 4lb.relabels_(["d","e"])
 5
 6
 7## contract all
 8X = cytnx.Contract(cytnx.Contract(A,la),cytnx.Contract(B,lb))
 9lb_l = lb.relabel("e","a")
10X = cytnx.Contract(lb_l,X)
11
12
13## X =
14#           (0)  (1)
15#            |    |
16#  (d) --lb-A-la-B-lb-- (e)
17#
18# X.print_diagram()
19Xt = X.clone()
20
21
22## calculate norm and energy for this step
23# Note that X,Xt contract will result a rank-0 tensor, which can use item() toget element
24XNorm = cytnx.Contract(X,Xt).item()
25XH = cytnx.Contract(X,H)
26XH.relabels_(["d","e","0","1"])
27
28
29XHX = cytnx.Contract(Xt,XH).item() ## rank-0
30E = XHX/XNorm
31
32# print(E)
33## check if converged.
34if(np.abs(E-Elast) < CvgCrit):
35    print("[Converged!]")
36    break
37print("Step: %d Enr: %5.8f"%(i,Elast))
38Elast = E

in the next step we perform the two-sites imaginary time evolution, using the operator (or “gate”) eH we defined above:

../_images/itebd_envolve.png

we also performed SVD for the XeH here, this put the MPS into mixed canonical form and have a Schimit decomposition of the whole state where the singular values are simply the Schimit coefficients. The Svd_truncate is called such that the intermediate bonds with label (-6) and (-7) are properly truncate to the maximum virtual bond dimension chi.

  • In Python

1## Time evolution the MPS
2XeH = cytnx.Contract(X,eH)
3XeH.permute_(["d","2","3","e"])
4XeH.set_rowrank_(2)
5la,A,B = cytnx.linalg.Svd_truncate(XeH,chi)
6la.normalize_()

Note that we directly store the SVD results into A, B and la, this can be seen by comparing to our original MPS configuration:

../_images/itebd_what.png

to recover to orignial form, we put \(\lambda_B^{-1} \lambda_B\) on both ends, which abosorb two \(\lambda_B^{-1}\):

../_images/itebd_recover.png

Now we have the envolved \(\Gamma_A\), \(\Gamma_B\) and \(\lambda_A\). Using the translation symmetry, we shift the whole chain to left by just exchange the \(Gamma\) and \(\lambda\) pair and arrived at the new MPS for next iteration to update B-A sites using \(U_b\).

../_images/itebd_translation.png

  • In Python

1lb_inv = 1./lb
2lb_inv.relabels_(["e","d"])
3A = cytnx.Contract(lb_inv,A)
4B = cytnx.Contract(B,lb_inv)
5# translation symmetry, exchange A and B site
6A,B = B,A
7la,lb = lb,la

Let’s put everything together in a loop for iteration:

  • In Python

 1for i in range(10000):
 2
 3    A.relabels_(["a","0","b"])
 4    B.relabels_(["c","1","d"])
 5    la.relabels_(["b","c"])
 6    lb.relabels_(["d","e"])
 7
 8    ## contract all
 9    X = cytnx.Contract(cytnx.Contract(A,la),cytnx.Contract(B,lb))
10    lb_l = lb.relabel("e","a")
11    X = cytnx.Contract(lb_l,X)
12    Xt = X.clone()
13
14    ## calculate norm and energy for this step
15    # Note that X,Xt contract will result a rank-0 tensor, which can use item() toget element
16    XNorm = cytnx.Contract(X,Xt).item()
17    XH = cytnx.Contract(X,H)
18    XH.relabels_(["d","e","0","1"])
19
20    XHX = cytnx.Contract(Xt,XH).item() ## rank-0
21    E = XHX/XNorm
22
23    # print(E)
24    ## check if converged.
25    if(np.abs(E-Elast) < CvgCrit):
26        print("[Converged!]")
27        break
28    print("Step: %d Enr: %5.8f"%(i,Elast))
29    Elast = E
30
31
32    ## Time evolution the MPS
33    XeH = cytnx.Contract(X,eH)
34    XeH.permute_(["d","2","3","e"])
35
36    ## Do Svd + truncate
37    XeH.set_rowrank_(2)
38    la,A,B = cytnx.linalg.Svd_truncate(XeH,chi)
39    la.normalize_()
40
41    lb_inv = 1./lb
42    lb_inv.relabels_(["e","d"])
43
44    A = cytnx.Contract(lb_inv,A)
45    B = cytnx.Contract(B,lb_inv)
46
47    # translation symmetry, exchange A and B site
48    A,B = B,A
49    la,lb = lb,la
[Vid07] (1,2)

G. Vidal. Classical simulation of infinite-size quantum lattice systems in one spatial dimension. Phys. Rev. Lett., 98:070201, Feb 2007. URL: https://link.aps.org/doi/10.1103/PhysRevLett.98.070201, doi:10.1103/PhysRevLett.98.070201.