# Different output each time I run a function

**URL:** <https://itensor.discourse.group/t/different-output-each-time-i-run-a-function/2364>\
**Category:** ITensor Julia Questions\
**Tags:** julia, mps\
**Created:** [April 29, 2025, 3:10pm UTC](https://itensor.discourse.group/t/different-output-each-time-i-run-a-function/2364 "2025-04-29T15:10:59Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![phx](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/phx/32/673_2.png) [@phx](https://itensor.discourse.group/u/phx)\
**Post date:** [April 29, 2025, 3:10pm UTC](https://itensor.discourse.group/t/different-output-each-time-i-run-a-function/2364/1 "2025-04-29T15:10:59Z")

</div>

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:

```auto
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 = N`then 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

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

```

five times, one after the other.

 ![test](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/4dc6d76d917de1224772b38a1bedc65b7dafad88.png)

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [April 29, 2025, 10:07pm UTC](https://itensor.discourse.group/t/different-output-each-time-i-run-a-function/2364/2 "2025-04-29T22:07:44Z")

</div>

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

```julia
rnd = 0 .* random_mps(...)

```

But still that’s not the right thing to do. You should use [`MPS` constructor](https://docs.itensor.org/ITensorMPS/stable/MPSandMPO.html#ITensorMPS.MPS-Tuple%7BType%7B%3C:Number%7D,%20Vector%7B%3C:Index%7D%7D)

```julia
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?

---

<div class="post-metadata">

**Author:** ![phx](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/phx/32/673_2.png) [@phx](https://itensor.discourse.group/u/phx)\
**Post date:** [April 30, 2025, 12:27am UTC](https://itensor.discourse.group/t/different-output-each-time-i-run-a-function/2364/3 "2025-04-30T00:27:42Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![system](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/d3072b13e047cd06df7f3981547c0917d940dfa4.png) [@system](https://itensor.discourse.group/u/system)\
**Post date:** [May 10, 2025, 12:27am UTC](https://itensor.discourse.group/t/different-output-each-time-i-run-a-function/2364/4 "2025-05-10T00:27:52Z")

</div>

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.