4. iDMRG

By : Hsu Ke, Kai-Hsin Wu

In the previous example we demostrated how the finite size DMRG is implemented, Here we present the infinite system size DMRG (iDMRG) scheme which consider the system with infinite boundary condition[PVM12], where the bondaries (enviroments) are grows the lattice by one or more sites each iteration, at the fixed point, a translationally invariant wavefunction is produced.

In another algorithm of infinite size simulation, iTEBD, a properly converged state requires careful scaling of the time-step to zero, a DMRG approach where an efficient local eigensolver is used to find a variationally optimal state ought to be much more robust and efficient. There are also disadavatages for iDMRG, we recommend the interested reader to refer to the paper of iDMRG or benchmarks made in the relevant papers.

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

\[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, refer to the previous examples.

The MPO is constucted as following:

\[\begin{split} M_j = \begin{bmatrix} I & S^z_j& H_x S^x_j \\ 0 & 0 & S^z_j \\ 0 & 0 & I \end{bmatrix}\end{split}\]

with the left and right boundary:

\[ \begin{align}\begin{aligned}\begin{split}L = \begin{bmatrix} 1\\ 0\\ 0 \end{bmatrix}\end{split}\\\begin{split}R = \begin{bmatrix} 0\\ 0\\ 1 \end{bmatrix}\end{split}\end{aligned}\end{align} \]

4.1. Initialization

The initailzation of MPO is much the same as we did in the previous DMRG example:

  • In Python:

 1J  = 1.0
 2Hx = 1.0
 3
 4d = 2
 5sx = cytnx.physics.pauli("x").real()
 6sz = cytnx.physics.pauli("z").real()
 7eye = cytnx.eye(d)
 8M = cytnx.zeros([3, 3, d, d])
 9M[0,0] = M[2,2] = eye
10M[0,1] = M[1,2] = sz
11M[0,2] = 2*Hx*sx
12M = cytnx.UniTensor(M,rowrank=0)
13
14L0 = cytnx.UniTensor(cytnx.zeros([3,1,1]),rowrank=0) #Left boundary
15R0 = cytnx.UniTensor(cytnx.zeros([3,1,1]),rowrank=0) #Right boundary
16L0.get_block_()[0,0,0] = 1.; R0.get_block_()[2,0,0] = 1.;
17
18L = L0
19R = R0

4.2. Update procedure

The first step of iDMRG is to obtain a initial two-site state, which is done by solving the eigenvalue problem of two-site hamitonian:

../_images/idmrg_n0.png

Let’s implement the function solving eigenvalue problem using in-built Lanczos method and the needed linear operation class:

  • In Python:

 1class Projector(cytnx.LinOp):
 2
 3
 4    def __init__(self,L,M1,M2,R,psi_dim,psi_dtype,psi_device):
 5        cytnx.LinOp.__init__(self,"mv",psi_dim,psi_dtype,psi_device)
 6
 7        self.anet = cytnx.Network("projector.net")
 8        self.anet.PutUniTensor("M2",M2)
 9        self.anet.PutUniTensors(["L","M1","R"],[L,M1,R])
10        self.psi_shape = [L.shape()[1],M1.shape()[2],M2.shape()[2],R.shape()[1]]
11
12    def matvec(self,psi):
13
14        psi_p = cytnx.UniTensor(psi.clone(),rowrank=0)  ## clone here
15        psi_p.reshape_(*self.psi_shape)
16
17        self.anet.PutUniTensor("psi",psi_p)
18        H_psi = self.anet.Launch().get_block_() # get_block_ without copy
19
20        H_psi.flatten_()
21        return H_psi
22
23
24
25def eig_Lanczos(psivec, functArgs, Cvgcrit=1.0e-15,maxit=100000):
26    """ Lanczos method for finding smallest algebraic eigenvector of linear \
27    operator defined as a function"""
28    #print(eig_Lanczos)
29
30    Hop = Projector(*functArgs,psivec.shape()[0],psivec.dtype(),psivec.device())
31    gs_energy ,psivec = cytnx.linalg.Lanczos(Hop,Tin=psivec,method="Gnd",CvgCrit=Cvgcrit,Maxiter=maxit)
32
33    return psivec, gs_energy.item()

Now do the optimization and SVD task:

  • In Python:

1psi = cytnx.UniTensor(cytnx.random.normal([1,d,d,1],1,2),rowrank=2)
2shp = psi.shape()
3psi_T = psi.get_block_(); psi_T.flatten_() ## flatten to 1d
4psi_T, Entemp = eig_Lanczos(psi_T, (L,M,M,R), maxit=maxit);
5psi_T.reshape_(*shp)
6psi = cytnx.UniTensor(psi_T,rowrank=2)
7
8s0,A,B = cytnx.linalg.Svd_truncate(psi,min(chi,d)) ## Svd
9s0/=s0.get_block_().Norm().item() ## normalize

Note

we are using a unit cell of two sites, however the unit cell can be any size, including a single site.

we performed SVD and use the left and right basis to update the environment for effective hamitonian, these procedure and network will always be the same in the future interations

../_images/idmrg_envupdate.png

  • In Python:

1anet = cytnx.Network("L_AMAH.net")
2anet.PutUniTensors(["L","A","A_Conj","M"],[L,A,A.Conj(),M]);
3L = anet.Launch()
4anet = cytnx.Network("R_AMAH.net")
5anet.PutUniTensors(["R","B","B_Conj","M"],[R,B,B.Conj(),M]);
6R = anet.Launch()

we then solve the eigenvalue problem again and do SVD for the new effective hamitonian, note that we initialized a new random trial state here.

../_images/idmrg_n1.png

  • In Python:

1## Construct n = 1
2psi = cytnx.UniTensor(cytnx.random.normal([d,d,d,d],0,2),rowrank=2)
3shp = psi.shape()
4psi_T = psi.get_block_(); psi_T.flatten_() ## flatten to 1d
5psi_T, Entemp = eig_Lanczos(psi_T, (L,M,M,R), maxit=maxit);
6psi_T.reshape_(*shp)
7psi = cytnx.UniTensor(psi_T,rowrank=2)
8s1,A,B = cytnx.linalg.Svd_truncate(psi,min(chi,d*d))
9s1/=s1.get_block_().Norm().item()

followed by another environment update:

  • In Python:

1# absorb A[1], B[1] to left & right enviroment.
2anet = cytnx.Network("L_AMAH.net")
3anet.PutUniTensors(["L","A","A_Conj","M"],[L,A,A.Conj(),M]);
4L = anet.Launch()
5anet = cytnx.Network("R_AMAH.net")
6anet.PutUniTensors(["R","B","B_Conj","M"],[R,B,B.Conj(),M]);
7R = anet.Launch()

The next few steps involve “rotate” the center of our state to the left and right:

../_images/idmrg_rotate.png

which is done by the straightforward contraction and re-SVD:

  • In Python:

1A.set_rowrank_(1)
2# set the label '_aux_R' to another to avoid Svd error
3sR,_,A = cytnx.linalg.Svd(cytnx.Contract(A, s1).relabel_(2, "-2"))
4B.set_rowrank_(2)
5# set the label '_aux_L' to another to avoid Svd error
6sL,B,_ = cytnx.linalg.Svd(cytnx.Contract(s1, B).relabel_(0, "-1"))
7A,B = B,A

Note that we have discarded U for the rotate left case and Vt for the right, and now in this scheme we construct our trial state for the n=2 iteration as follows, which looks like a inserted two new sites between our environment:

../_images/idmrg_trial.png

and optimize this state:

../_images/idmrg_n2.png

But where does that trial form come from? To gain more insight about it, we can instead consider the iTEBD approach here, which is done by perform an evolution “gate” on the two sites:

../_images/idmrg_gate.png

Suppose we didn’t even know what \(s_n\) looks like here, it’s just a \(s_n\), then we can simply follows the procedures mentioned in the iTEBD example, we have the following Vidal’s form for unit cell of two sites at n-th step:

../_images/idmrg_vidal.png

this immediately suggests the above trial form if we try to rotate the center to the left and right and use translation symmetry.

The construction of trial state and optimization is done as follows:

  • In Python:

 1sR.relabel_(0,"1")
 2sL.relabel_(1,"0")
 3s0 = 1./s0
 4s0.relabels_(["0","1"])
 5s2 = cytnx.Contract(cytnx.Contract(sL,s0),sR)
 6
 7s2.relabels_(["-10","-11"])
 8A.relabel_(2,"-10")
 9B.relabel_(0,"-11")
10psi = cytnx.Contract(cytnx.Contract(A,s2),B)
11
12## optimize wave function:
13#  again use Lanczos to get the (local) ground state, the projector is in shape
14shp = psi.shape()
15psi_T = psi.get_block_(); psi_T.flatten_() ## flatten to 1d
16psi_T, Entemp = eig_Lanczos(psi_T, (L,M,M,R), maxit=maxit);
17psi_T.reshape_(*shp)
18psi = cytnx.UniTensor(psi_T,rowrank=2)
19s2,A,B = cytnx.linalg.Svd_truncate(psi,min(chi,psi.shape()[0]*psi.shape()[1]))
20s2/=s2.get_block_().Norm().item()

then we check the convergence by comparing the singular values to the one we obtained in the previous iteration:

  • In Python:

1if(s2.get_block_().shape()[0] != s1.get_block_().shape()[0]):
2    ss = 0
3    print("step:%d, increasing bond dim!! dim: %d/%d"%(i,s1.get_block_().shape()[0],chi))
4else:
5    ss = abs(cytnx.linalg.Dot(s2.get_block_(),s1.get_block_()).item())
6    print("step:%d, diff:%11.11f"%(i,1-ss))
7if(1-ss<1.0e-10):
8    print("[converge!!]")
9    break;

also rememeber to update the environment using the SVD result.

  • In Python:

 1anet = cytnx.Network("L_AMAH.net")
 2anet.PutUniTensors(["L","A","A_Conj","M"],[L,A,A.Conj(),M]);
 3L = anet.Launch()
 4
 5anet = cytnx.Network("R_AMAH.net")
 6anet.PutUniTensors(["R","B","B_Conj","M"],[R,B,B.Conj(),M]);
 7R = anet.Launch()
 8
 9s0 = s1
10s1 = s2

After reaching the fixed point, let’s consider a local measurement of energy for the final state:

../_images/idmrg_measure.png

  • In Python:

1H = J*cytnx.linalg.Kron(sz,sz) + Hx*(cytnx.linalg.Kron(sx,eye) + cytnx.linalg.Kron(eye,sx))
2H = cytnx.UniTensor(H.reshape(2,2,2,2),rowrank=2)
3
4# use the converged state to get the local energy:
5anet = cytnx.Network("Measure.net")
6anet.PutUniTensors(["psi","psi_conj","M"],[psi,psi,H])
7E = anet.Launch().item()
8print("ground state E",E)
[PVM12]

Ho N. Phien, Guifré Vidal, and Ian P. McCulloch. Infinite boundary conditions for matrix product state calculations. Phys. Rev. B, 86:245107, Dec 2012. URL: https://link.aps.org/doi/10.1103/PhysRevB.86.245107, doi:10.1103/PhysRevB.86.245107.