Product MPOs have bond dimension greater than 1 with OpSum

[mtfishman]

@emstoudenmire, I think this issue has come up before (but maybe it was in the C++ version). Here is a minimal example:

using ITensors
s = siteinds("S=1/2", 3)
M = MPO(AutoMPO() + ("Sz", 1), s)
@show inds(M[2])

returns:

inds(M[2]) = IndexSet{4} (dim=3|id=714|"Link,l=1") (dim=3|id=687|"Link,l=2") (dim=2|id=144|"S=1/2,Site,n=2")' (dim=2|id=144|"S=1/2,Site,n=2")

It seems like any single term AutoMPO I input returns an MPO with bond dimension 3. I assume that is because it is making an MPO of the form that could accommodate other MPO terms, and doesn't specialize to the case of just a single term.

The reason this came up is that I want to implement a generic MPO constructor MPO(s, "Id") which in general returns an MPO that is a product of specified operators. This is currently available for non-QN MPOs, but as I was generalizing to QN MPOs I realized it would be a lot easier to let AutoMPO generate the MPO for me (since it already deals with all of the subtleties of dealing with blocks, flux, etc.). However, I noticed it doesn't output the MPO with the optimal bond dimension as shown above.

In theory we could truncate the MPO output by AutoMPO, but naively truncating an MPO with our truncate! function doesn't work in general, since the norm can diverge (like for the case of MPO(s, "Id"), the norm of that operator grows as 2^N). As a separate issue, we should remember to implement the MPO truncation algorithm from the paper https://arxiv.org/abs/1909.06341.

[mtfishman]

I know Miles is aware, but just for general knowledge and as a reminder, this is what happens if you try to naively truncate a Hamiltonian-like MPO:

julia> s = siteinds("S=1/2", 200);

julia> M = MPO(AutoMPO() + ("Id", 1), s);

julia> truncate!(M; cutoff = 1e-15);

julia> @show M[1];
M[1] = ITensor ord=3
Dim 1: (dim=2|id=375|"S=1/2,Site,n=1")'
Dim 2: (dim=2|id=375|"S=1/2,Site,n=1")
Dim 3: (dim=1|id=662|"Link,l=1")
NDTensors.Dense{Float64,Array{Float64,1}}
 2×2×1
[:, :, 1] =
 -8.96364335596578e29   0.0
  0.0                  -8.963643355965783e29

It does indeed truncate properly, but because it is done naively by treating the MPO as a state, the norm of one of the tensors explodes.

[mtfishman]

A temporary workaround to truncate these types of MPOs without a native AutoMPO fix is the following:

s = siteinds("S=1/2", 200)
M = MPO(AutoMPO() + ("Id", 1), s)
lognormM = lognorm(M)
M ./= exp(lognormM / N)
truncate!(M; cutoff = 1e-15)
M .*= exp(lognormM / N)

so the MPO is normalized using lognorm, which is numerically well behaved even for MPOs with diverging norms like the identity, and then truncated.

Ideally, we would have this as a setting in truncate!, where the user can input what type of truncation they want to perform depending on the type of MPO that is input. The default could be that it is treated as a density matrix-like MPO, where either a standard state-like truncation is used or a more sophisticated algorithm like DMT (https://arxiv.org/abs/1707.01506) could be used. Then, if it is a local Hamiltonian-like MPO, it could normalize the MPO using the norm-per-site as above, and in a more sophisticated version use the truncation from https://arxiv.org/abs/1909.06341.

[emstoudenmire]

I will think about this some more, but at a first pass over the AutoMPO code and logic, the case of bond dimension =1 is a rather special case, to the point where if AutoMPO was to handle it, it might be better to have a separate backend code inside AutoMPO just for this case. One way to see why it's special is that normally there are starting and ending identity strings which can be built from any MPO tensors, whereas for the product operator case, the starting and ending strings are merged together, and operators in them which are not identity are "statically" known ahead of time and merged in too (by merged I mean not given their own row or column in the MPO matrix).

[emstoudenmire]

If we do just make a special-case code, it would be a good idea to put it inside AutoMPO though so that way it will give the capability to both AutoMPO and the MPO constructors which draw upon it.

[mtfishman]

That makes sense that it would require a special code path within AutoMPO. Agreed that ideally this would be within AutoMPO.

Something to consider for this case of a product MPO is that we could just have it output an MPO without link indices, in which case the implementation would be pretty trivial. I have been meaning to make sure all of the MPS/MPO code "just works" if there are missing link indices, and this would be a good use case for that.

Something else I wanted to bring up is that this functionality bumps into the limitation that AutoMPO doesn't handle the case where the MPO has a nonzero flux, i.e. we would want to make sure this works for things like MPO(s, "S+").