Modified TEBD with QR decomposition

ITensor Julia Questions
, ,
1

Hello everyone!

This paper [ Fast time evolution of matrix product states using the QR decomposition | Phys. Rev. B] develops a modified TEBD or Trotter gates time-evolving method in which SVD is replaced by QR+LQ decomposition. This modification accelerate the time evolution siginificantly campared to normal ones. Its main modification is presented as below.

image
image675×583 42.4 KB

However, this modified TEBD is applied to infinite MPS. Following the process posted before, I programed my version of this modified TEBD for center-orthogonalized finite MPS , in which the SVD for two updating tensors connected by the orthogonal center bond is replaced by QR+LQ decomposition. But the time-evolution results of the modified one and orignal one are particularly different.

After a lot efforts, I still can not figure it out. Can anyone give some suggestions? Should I transform the center-orthogonalized MPS to right-isometric MPS? How can ITensor construct a finite MPS in isometric form?

Here is my code

using ITensors
using ITensorMPS
using LinearAlgebra

### Trotter gates
L=6
tstep=0.1
gates=ITensor[]
sites=siteinds("S=1/2",L)
for i=1:L-1
    s1=sites[i]
    s2=sites[i+1]
    ops = -1/2*op("S+",s1)*op("S-",s2)
    ops += -1/2*op("S-",s1)*op("S+",s2)
    ops += op("Sz",s1)*op("Sz",s2)
    tevol=exp(-im*tstep*ops)
    push!(gates,tevol)
end


### Hamiltonian
hal=OpSum()
for i=1:L-1
    hal += -1/2,"S+",i,"S-",i+1
    hal += -1/2,"S-",i,"S+",i+1
    hal += "Sz",i,"Sz",i+1
end
Hal=toMPO(hal,sites)

### updating edge tensors with svd in open boudary condition
function tevol_edge(tensor_1::ITensor,tensor_2::ITensor,gate_12::ITensor,site_ind::Int64; kwargs...)
    s1s2_t1=tensor_1*tensor_2*gate_12
    link_ind=commonind(tensor_1,tensor_2)
    if site_ind==1
        outer_ind=setdiff(uniqueinds(tensor_2,tensor_1),[sites[site_ind+1]])[1]
        U1,S1,V1=svd(s1s2_t1,prime(sites[site_ind]);kwargs...)
        new_tags="Link,l="*string(site_ind)
        # U1=permute(U1,prime(sites[site_ind]),commonind(U1,S1))
        new_tensor_1=replaceind(U1,commonind(U1,S1),settags(commonind(U1,S1),new_tags))
        new_tensor_1=noprime(new_tensor_1,prime(sites[site_ind]))
        new_tensor_1=permute(new_tensor_1,sites[site_ind],settags(commonind(U1,S1),new_tags))
        # @show inds(new_tensor_1)
        new_tensor_2=replaceind(S1*V1,commonind(U1,S1),settags(commonind(U1,S1),new_tags))
        # @show commonind(new_tensor_2,new_tensor_1)
        new_tensor_2=noprime(new_tensor_2,prime(sites[site_ind+1]))
        # @show new_tensor_2
        new_tensor_2=permute(new_tensor_2,commonind(new_tensor_2,new_tensor_1),sites[site_ind+1],outer_ind)
    else
        outer_ind=setdiff(uniqueinds(tensor_1,tensor_2),[sites[site_ind]])[1]
        U1,S1,V1=svd(s1s2_t1,prime(sites[site_ind]),outer_ind;kwargs...)
        new_tags="Link,l="*string(site_ind)
        # U1=permute(U1,prime(sites[site_ind]),commonind(U1,S1))
        new_tensor_1=replaceind(U1,commonind(U1,S1),settags(commonind(U1,S1),new_tags))
        new_tensor_1=noprime(new_tensor_1,prime(sites[site_ind]))
        new_tensor_1=permute(new_tensor_1,outer_ind,sites[site_ind],settags(commonind(U1,S1),new_tags))
        # @show inds(new_tensor_1)
        new_tensor_2=replaceind(S1*V1,commonind(U1,S1),settags(commonind(U1,S1),new_tags))
        # @show commonind(new_tensor_2,new_tensor_1)
        new_tensor_2=noprime(new_tensor_2,prime(sites[site_ind+1]))
        # @show new_tensor_2
        new_tensor_2=permute(new_tensor_2,commonind(new_tensor_2,new_tensor_1),sites[site_ind+1])

    end

    return (new_tensor_1,new_tensor_2)
end

psi_rand=randomMPS(sites,linkdims=4)
orthogonalize!(psi_rand,1)
psi_rand_c=copy(psi_rand)

### updating site 1 and site 2 by SVD
t1=tevol_edge(psi_rand[1],psi_rand[2],gates[1],1;cutoff=0,maxdim=20)
psi_rand_c[1]=t1[1]
psi_rand_c[2]=t1[2]

### updating site 2 and site 3 by QR+LQ
B_m=psi_rand_c[2]
B_n=psi_rand_c[3]

theta_1=(B_m*B_n)*gates[2]*dag(prime(B_n,sites[3]))
Q_m,R=qr(theta_1,commonind(psi_rand_c[1],B_m),prime(sites[2]))
theta_2=(B_m*B_n)*gates[2]*dag(Q_m)
L,Q_n=lq(theta_2,commonind(Q_m,R))

#### update site 3
psi_rand_c[3]=replaceind(noprime(L*Q_n,prime(sites[3])),uniqueind(L,Q_n),settags(uniqueind(L,Q_n),"Link,l=2"))
psi_rand_c[3]=permute(psi_rand_c[3],settags(uniqueind(L,Q_n),"Link,l=2"),sites[3],commonind(psi_rand_c[3],psi_rand_c[4]))
gamma_m,B_m_1=lq(B_m,commonind(psi_rand_c[1],B_m))

#### update site 2
B_m_tilde=noprime((B_m_1*B_n)*gates[2]*dag(Q_n),prime(sites[2]))
B_m_tilde=replaceind(B_m_tilde,commonind(B_m_1,gamma_m),uniqueind(gamma_m,B_m_tilde))
B_m_tilde=replaceind(B_m_tilde,commonind(L,Q_n),settags(uniqueind(L,Q_n),"Link,l=2"))
psi_rand_c[2]=B_m_tilde

normalize!(psi_rand_c)
expect(psi_rand_c,"Sz")


### orignal TEBD with SVD
psi_rand_c1=copy(psi_rand)
psi_rand_c1[1]=t1[1]
psi_rand_c1[2]=t1[2]

theta_1_1=(psi_rand_c1[2]*psi_rand_c1[3])*gates[2]
new_B_m,S,new_B_n=svd(theta_1_1,(commonind(psi_rand_c1[1],psi_rand_c1[2]),prime(sites[2])))
new_B_n_1=S*new_B_n
new_B_n_1=replaceind(new_B_n_1,commonind(new_B_m,S),settags(commonind(new_B_m,S),"link,l=2"))
new_B_n_1=replaceind(new_B_n_1,prime(sites[3]),sites[3])
new_B_n_1=permute(new_B_n_1,settags(commonind(new_B_m,S),"link,l=2"),sites[3],commonind(new_B_n_1,psi_rand_c1[3]))
psi_rand_c1[3]=new_B_n_1
new_B_m=replaceind(new_B_m,prime(sites[2]),sites[2])
new_B_m=replaceind(new_B_m,commonind(new_B_m,S),settags(commonind(new_B_m,S),"link,l=2"))
psi_rand_c1[2]=new_B_m

normalize!(psi_rand_c1)
expect(psi_rand_c1,"Sz")

The expectation value of $$Sz$$ is particular different.

2

Hi Wang,
My main suggestion would be to try to work out the steps of the algorithm for a finite MPS. I have personally found that the finite MPS case can be an easier setting to understand and debug things versus the infinite setting which can be more technical.

As far as how ITensor represents isometric MPS, I’m not sure I totally understand the question there, but an isometric gauge (which I usually call an “orthogonal gauge”) of an MPS is a property of each tensor one by one from the left and right ends of the MPS, and we just store tensors with those properties in an MPS object in ITensor when we have brought an MPS into such a gauge. We have functions such as psi = orthogonalize(psi,j) which brings any finite MPS into an orthogonal gauge with the orthogonality “center” located at site j.