[Julia] Julia DMRG gives inconsistent result for weak interactions

+1 vote
edited

Hello,

It seems that the DMRG in Julia gives inconsistent result for weak interactions. The sample code is below.

ITensors.space(::SiteType"HardCore") = 2

ITensors.state(::StateName"0", ::SiteType"HardCore") = [1.0, 0.0]
ITensors.state(::StateName"1", ::SiteType"HardCore") = [0.0, 1.0]

function ITensors.op!(Op::ITensor,
::OpName"N",
::SiteType"HardCore",
s::Index)
Op[s'=>2,s=>2] = 1
end

function ITensors.op!(Op::ITensor,
::SiteType"HardCore",
s::Index)
Op[s'=>1,s=>2] = 1
end

function ITensors.op!(Op::ITensor,
::OpName"A",
::SiteType"HardCore",
s::Index)
Op[s'=>2,s=>1] = 1
end

t = 1.0
V1 = 1.0e-3
V2 = 2e-5

N = 20
sites = siteinds("HardCore",N)

ampo = OpSum()

# NN
for j = 1:N-1
ampo += V1,"N",j,"N",j+1
end

# NNN
for j = 1:N-2
ampo += V2,"N",j,"N",j+2
end

H = MPO(ampo,sites)

psi0 = productMPS(sites, n -> isodd(n) ? "0" : "1")
println(inner(psi0,H,psi0))

# Create an initial random matrix product state
psi0 = randomMPS(sites)

sweeps = Sweeps(100)
setmaxdim!(sweeps, 20,20,40,40,60,60,80)
setcutoff!(sweeps, 1E-11)
#noise!(sweeps, 1E-4,1E-6,1E-10,0)

N_observer = DMRGObserver(["N"],sites,energy_tol=1E-12)

energy, psi = dmrg(H,psi0, sweeps, observer=N_observer)
println("Final energy = \$energy")


The result from the code is -12.378711861320. But exact diagonalization gives -12.378713430075452. There is a 1.6E-6 discrepancy. Although it is small, the DMRG in C++ version gives -12.378713429701, only 3.7e-10 discrepancy with ED. The sample code in C++ is below.

#include "itensor/all.h"

using namespace std;
using namespace itensor;

int
main(int argc, char* argv[])
{
auto Ns = 20;
auto t = 1.0;
auto V1 = 1E-3;
auto V2 = 2E-5;

auto sites = Boson(Ns,{"ConserveQNs=",true});

auto ampo = AutoMPO(sites);

for(auto j : range1(Ns-1))
{
ampo += V1,"N",j,"N",j+1;
}

for(auto j : range1(Ns-2))
{
ampo += V2,"N",j,"N",j+2;
}

auto H = toMPO(ampo);

//Perform 20 sweeps of DMRG
auto sweeps = Sweeps(20);
//Specify max number of states kept each sweep
sweeps.cutoff() = 1E-11;
sweeps.maxdim() = 20, 20, 40, 40, 60;

auto state = InitState(sites);
for(auto j : range1(Ns))
{
state.set(j, j%2 == 1 ? "Emp" : "1");
}
auto psi0 = MPS(state);

auto [energy0, psi0f] = dmrg(H, psi0, sweeps, {"Quiet=", true, "EnergyErrgoal=",1E-12,"EntropyErrgoal=",1E-12});
cout << energy0 << endl;

return 0;
}


What causes this difference in Julia? Thanks.

Jin

commented by (70.1k points)
It’s hard to know what caused the small difference, but I note that both DMRG calculations are approximate and have different levels of approximation (in this case, different final maxdim). Could you please try some DMRG runs with even larger final maxdim, more sweeps, and a smaller cutoff to verify that the results become closer to ED?
commented by (1.2k points)
Hi Miles, thanks for the quick comment.
When I set cutoff = 1E-15, final maxdim=100 in Julia, I got the following results
After sweep 1 energy=-12.239362658662 maxlinkdim=4 maxerr=0.00E+00 time=0.017
After sweep 2 energy=-12.373944285627 maxlinkdim=16 maxerr=0.00E+00 time=0.054
After sweep 3 energy=-12.378703148154 maxlinkdim=40 maxerr=9.71E-11 time=0.136
After sweep 4 energy=-12.378711857431 maxlinkdim=40 maxerr=1.08E-10 time=0.185
After sweep 5 energy=-12.378711861708 maxlinkdim=60 maxerr=1.08E-13 time=0.290
After sweep 6 energy=-12.378711861708 maxlinkdim=60 maxerr=9.86E-14 time=0.378
After sweep 7 energy=-12.378711861709 maxlinkdim=80 maxerr=5.16E-15 time=0.394
After sweep 8 energy=-12.378711861709 maxlinkdim=80 maxerr=5.11E-15 time=0.411
After sweep 9 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.414
After sweep 10 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.673
After sweep 11 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.538
After sweep 12 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.501
After sweep 13 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.426
After sweep 14 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.437
After sweep 15 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.424
After sweep 16 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.494
After sweep 17 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.505
After sweep 18 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.752
After sweep 19 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.569
After sweep 20 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.500
After sweep 21 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.420
After sweep 22 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.460
After sweep 23 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.451
After sweep 24 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.425
After sweep 25 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.418
After sweep 26 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.426
After sweep 27 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.505
After sweep 28 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.474
After sweep 29 energy=-12.378711861709 maxlinkdim=84 maxerr=9.90E-16 time=0.632
Energy difference less than 2.0e-15, stopping DMRG
The result still has 1.6E-6 discrepancy.

For C++ version with the same parameter, the final step reads
vN Entropy at center bond b=10 = 0.757996415828
Eigs at center bond b=10: 0.7680 0.1077 0.1077 0.0151
Largest link dim during sweep 16/100 was 84
Largest truncation error: 7.6029e-16
Energy after sweep 16/100 is -12.378713430075
Sweep 16/100 CPU time = 0.148s (Wall time = 0.154s)
Energy error dE = 0.000E+00 after 16 sweeps;
Entropy error dS = 9.992E-15 after 16 sweeps;
Energy (Entropy) error goal met (dE = 0.000E+00 < 2.000E-15 and dS = 9.992E-15 < 1.000E-12); returning after 16 sweeps.

The result from DMRG in C++ version is exactly the same as ED.
commented by (1.2k points)
If I make an initial state |0101...01> with 20 sites, the expectation value of H should be <H> = 9V_2 = 0.00018. The C++ code gives exact 0.00018. But the Julia code gives 0.00017992806074756208. The difference is much larger than the machine precision. It seems there's something wrong in the MPO of the Hamiltonian.
commented by (70.1k points)
Although it's possible, and we don't know the exact reason for this small difference yet, I would bet it is not due to the MPO. Most likely, it is due to some subtle convergence parameter settings within the "eigensolver" subroutine within the Julia DMRG code. This routine has various parameters which adjust how much of an accuracy it strives for, and it may be that we need to adjust these for very high-accuracy calculations like yours.

Thanks for noticing this issue - I will file a bug report so that we make sure to look into it.

Is it urgent or blocking any work you need to do?

Thanks,
Miles
commented by (70.1k points)
commented by (1.2k points)
Thanks for the comment. But in the last comment I didn't do DMRG, I just constructed MPO and measured the energy with respect to the state |0101...01>:  <0101...01| H |0101...01> = 0.00017992806074756208. The exact value should be 9V_2 = 0.00018. Isn't this due to incorrect MPO?

Best,
Jin
commented by (70.1k points)
Thanks for clarifying this - I missed this on my first read. Yes, that definitely looks like an inaccuracy in the MPO, so we will look into it and I've update the bug report that I filed to note that it could be an issue with the AutoMPO / OpSum system.

+1 vote

Thanks for catching this issue. I have just created a fix for the OpSum / AutoMPO system that makes the Hamiltonian more accurate in this case.

https://github.com/ITensor/ITensors.jl/pull/726

Please let me know if the other issue persists, where even after many DMRG sweeps the result does not match ED closely enough.

Miles

commented by (1.2k points)
Hi Miles, thanks for doing this. When using a smaller cutoff 1E-15 for MPO, the energy for the state |0101...01> is 0.00018000000000000007, very accurate. Then the DMRG result for the ground state energy is -12.378713430075, same as ED. So the DMRG has no problem. When I decrease V2 from 2E-5 to 2E-6, the issue happens again. The energy for the state |0101...01> is 1.7999928000607997e-5, different from the exact value 1.8E-5. The C++ version also gives 1.799992800e-5. So C++ also uses a cutoff 1E-15 to compress MPO, right?
commented by (1.2k points)
I thought the MPO representation is exact for short range interactions. Why do we do compression for MPO?
commented by (1.2k points)
I found C++ version uses cutoff = 1E-13, but is more accurate than Julia. Maybe the compression procedure has some difference.
commented by (70.1k points)
It's a good question about the compression. The reasons are a bit complicated, but basically you can imagine that the MPO starts out with matrices inside which are sparse and have a pattern that is easier for a human to write a code to make, but as a result are sometimes larger than they have to be. So then we use an essentially exact compression to find low-rank structure in these matrices, throwing away only singular values which are essentially zero. But it's that question about what "essentially zero" means which is subtle.

Good point about the C++ using 1E-13. Not sure why the result is different, but I think it may have to do with subtle differences in the function called truncate! in the ITensors.NDTensors module and the related function in C++.
commented by (1.2k points)
Thanks, Miles.