+1 vote
asked by (280 points)

Hi there,

I am trying to use the fitApplyMPO method, however it seems to me that I only find a local minimum. My start state is a product state and after an arbitrary number of sweeps, I always get back a product state. The bond dimension of my IQMPS does not increase.

Then I saw that your fitApplyMPO also allows for a noise term, that is added when creating the svdBond. I gave it a try and it always crashes with a segmentation fault.

I was trying to understand what is going on and I am not so sure that a noise term makes sense at all here.

To my understanding, you need the "environment" to have a paired index, which is not the case anymore. In the fitApply concept, the environmental tensors have three bond indices: one from the MPS you are optimizing, one from the MPO, and one from the MPS you are applying the MPO to. The two MPS links can be completely different, which then causes the fitApplyMPO method to crash.

Is my understanding right? Does anybody have an idea how to fix that?

Lars

1 Answer

0 votes
answered by (70.1k points)

Hi Lars,
That feature of fitApplyMPO may not be working. However, I recently wrote a slightly different version of fitApplyMPO that does support a noise term and I am pasting it below. Could you look through it and give it a try and let me know if it works for you (perhaps with minor modifications if there's a minor bug somewhere)?

If it works well for you I may replace what's in the library with it:

namespace itensor {

template<typename T>
MPSt<T>
fitApplyOneSite(MPSt<T> const& psi,
                MPOt<T> const& W,
                Sweeps const& sweeps,
                Args args = Args::global())
    {
    auto quiet = args.getBool("Quiet",true);
    auto N = psi.N();

    auto res = psi;

    auto LR = std::vector<T>(N+2);
    LR.at(N) = psi.A(N) * W.A(N) * dag(prime(res.A(N)));
    for(auto j = N-1; j > 1; --j)
        {
        LR.at(j) = LR.at(j+1)*psi.A(j)*W.A(j)*dag(prime(res.A(j)));
        }

    for(auto sw : range1(sweeps.nsweep()))
        {
        args.add("Sweep",sw);
        args.add("Cutoff",sweeps.cutoff(sw));
        args.add("Minm",sweeps.minm(sw));
        args.add("Maxm",sweeps.maxm(sw));
        auto noise = sweeps.noise(sw);
        for(int j = 1, ha = 1; ha <= 2; sweepnext1(j,ha,N))
            {
            if(not quiet) printfln("sw=%d j=%d ha=%d noise=%.3E",sw,j,ha,noise);

            //
            // Apply MPO
            //
            auto WA = psi.A(j);
            if(LR.at(j-1)) WA *= LR.at(j-1);
            WA *= W.A(j);
            if(LR.at(j+1)) WA *= LR.at(j+1);
            res.setA(j,noprime(WA));

            if(not (ha==2 && j==1))
                {
                res.Anc(j) /= norm(res.A(j));

                //
                // Form density matrix
                //
                auto d = ha==1 ? +1 : -1;
                auto li = linkInd(res,ha==1 ? j : j-1);
                auto rho = li ? primeExcept(res.A(j),li) : prime(res.A(j));
                rho *= dag(res.A(j));

                //
                // Form noise * density matrix perturbation
                //
                if(noise > 0.)
                    {
                    auto drho = psi.A(j);
                    if(LR.at(j-d)) drho *= LR.at(j-d);
                    drho *= W.A(j);
                    drho = dag(noprime(drho)) * drho;

                    rho += noise*drho;
                    }

                //
                // Diag density matrix
                //
                T U,D;
                auto spec = diagHermitian(rho,U,D,args);
                auto ci = commonIndex(U,D);
                if(not quiet) printfln("m=%d truncerr=%.5E",ci.m(),spec.truncerr());

                if(j+d <= N && j+d >= 1)
                    {
                    auto C = dag(U)*res.A(j);
                    res.setA(j,U);
                    res.setA(j+d,C*res.A(j+d));
                    }

                LR.at(j) = LR.at(j-d) ? LR.at(j-d)*psi.A(j) : psi.A(j);
                LR.at(j) *= W.A(j);
                LR.at(j) *= dag(prime(res.A(j)));
                }
            }
        }

    return res;
    }



} //namespace itensor

The above code is intended to go into a header file. You can put it into a .cc file too with using namespace itensor; above and just remove the namespace itensor { } part around it.

Best,
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.

Categories

...