Different output each time I run a function

ITensor Julia Questions
,
1

Hi everyone,
I’m struggling with the following problem: I have a function that manipulates an MPS and gives me a different result each time I execute it.

This is a minimal working example:

using ITensors
using ITensorMPS
using ITensorMPS: set_nsite!, position!, setleftlim!, setrightlim!
using KrylovKit: exponentiate

function solver(H, dt, state_0; kwargs...)
    return exponentiate(
        H,
        dt,
        state_0;
        tol=get(kwargs, :solver_tol, 1E-12),
        krylovdim=get(kwargs, :solver_krylovdim, 30),
        maxiter=get(kwargs, :solver_maxiter, 100),
        verbosity=get(kwargs, :solver_outputlevel, 0),
        eager=true,
    )
end

function incomplete_tdvp(N, dt, maxt, m=N; kwargs...)
    coups = ones(N - 1)
    freqs = [1; 0.5 * ones(N - 1)]

    sites = siteinds("S=1/2", N)

    state_t = MPS(sites, n -> n == 1 ? "Up" : "Dn")
    # Artificially increase state_t's bond dimensions
    rnd = 0 * random_mps(sites; linkdims=get(kwargs, :maxdim, 1) - 1)
    enlarged_state_t = add(state_t, rnd; alg="directsum")
    @assert isapprox(dot(state_t, enlarged_state_t), 1)
    state_t = enlarged_state_t

    h = OpSum()
    h += freqs[1], "Sz", 1
    for j in 2:N
        h += coups[j - 1], "S+", j - 1, "S-", j
        h += coups[j - 1], "S-", j - 1, "S+", j
        h += freqs[j], "Sz", j
    end
    PH = ProjMPO(MPO(h, sites))

    nsteps = floor(Int, maxt / dt)
    t = zero(dt):dt:maxt
    Sz = Vector{Float64}(undef, nsteps + 1)
    Sz[1] = expect(state_t, "Sz"; sites=1)

    for step in 1:nsteps
        orthogonalize!(state_t, 1)
        for site in 1:m  # if m != N it stops before the end of the MPS
            set_nsite!(PH, 1)
            position!(PH, state_t, site)

            # Forward evolution
            evolved_1site_tensor, info = solver(PH, -im * dt, state_t[site])
            info.converged == 0 && throw("exponentiate did not converge")

            if site < m  # Also backward evolution
                U, C = qr(
                    evolved_1site_tensor,
                    uniqueinds(evolved_1site_tensor, state_t[site + 1]);
                    tags="Link,n=$site",
                )
                state_t[site] = U
                setleftlim!(state_t, site)

                # Prepare the zero-site projection.
                set_nsite!(PH, 0)
                position!(PH, state_t, site + 1)

                C, info = solver(PH, im * dt, C)

                # Reunite the backwards-evolved C with the matrix on the next site.
                state_t[site + 1] *= C

                setrightlim!(state_t, site + 2)  # Now the orthocenter is on `site+1`
            else  # site == m
                # Nothing to do if the sweep is at the last site.
                state_t[site] = evolved_1site_tensor
            end
        end

        Sz[step + 1] = expect(state_t, "Sz"; sites=1)
    end
    return t, Sz
end

The incomplete_tdvp function basically performs a time evolution with a very crude implementation of the TDVP1 algorithm. I don’t really care here whether it’s correct or not from a physical point of view, as I’m using a more refined function in “real life” anyway.
Now, if the function gets called with m = Nthen it’s a complete TDVP left-to-right sweep and nothing bad happens.
However, when m < N the function behaves in a weird way, because each time I run it I get a different output.

I’m fine with a different output with respect to the m = N case, of course it’s going to be different (again, it doesn’t have to make sense physically) but I can’t see why it returns a different Sz anytime.

For example, here’s a plot generated by running using Plots and then

plot!(incomplete_tdvp(10, 0.1, 10, 4; cutoff=1e-10, maxdim=10))

five times, one after the other.

test
test600×400 30.6 KB

2

I think the issue is with how you are defining rnd, since scalar * MPS does not work how you think it does, you probably wanted to do

rnd = 0 .* random_mps(...)

But still that’s not the right thing to do. You should use MPS constructor

rnd = MPS(sites; linkdims=get(kwargs, :maxdim, 1) - 1)

By changing this line, this works. I have no idea on why it didn’t give you problems with m == N.

On a side note, I’m not sure if there is a way to avoid completely that part.

Explanation on * vs .*

when you do scalar * mps, the orthogonality center of the mps is rescaled by that factor, so in your case, only the first tensor of the random_mps was changing, the rest was random so it was different between calls. By doing .* instead, each tensor of mps gets multiplied by that scalar.
@mtfishman about 0*mps should a check with a warning/error be added?

2 Likes
3

You’re absolutely right, 0 * rnd calls some _apply_to_orthocenter! function that effectively leaves the sites of the MPS that are not the orthocenter untouched. In fact iszero.(0 * rnd) gives [1, 0, ..., 0] whereas iszero.(0 .* rnd) returns [1, 1, ..., 1].

When m == N maybe, since the sweep changes all the blocks of the MPS, these “fake zero” blocks get transformed into something that is effectively zero and does not affect the expectation value.

In any case, by using 0 .* rnd the output is now consistent. Thank you very much for spotting the issue!

1 Like