+1 vote
asked by (1.2k points)
edited by


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,
    Op[s'=>2,s=>2] = 1

function ITensors.op!(Op::ITensor,
    Op[s'=>1,s=>2] = 1

function ITensors.op!(Op::ITensor,
    Op[s'=>2,s=>1] = 1

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 += -t,"Adag",j,"A",j+1
    ampo += -t,"A",j,"Adag",j+1
    ampo += V1,"N",j,"N",j+1

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

H = MPO(ampo,sites)

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

# 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;

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 += -t,"Adag",j,"A",j+1;
        ampo += -t,"A",j,"Adag",j+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.


commented by (64.7k 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 (64.7k 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?

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?

commented by (64.7k 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 Answer

+1 vote
answered by (64.7k points)

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.


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


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 (64.7k 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.
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.