Product MPOs have bond dimension greater than 1 with OpSum

@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.