Can you use a noise term with fitApplyMPO

Question

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

Answer 1

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