[ITensors] [BUG] Inconsistent Behavior of mindim

[ryanlevy]

Description of bug

When using mindim!, the link dimension behavior at the "edge" is different depending on the value of cutoff!.
If there is a cutoff, some tensors get set to the mindim while others are not, and if there isn't a cutoff (cutoff!(sweeps,0)) the link dimensions are symmetric. I'm not sure which output is preferred but I suspect its the latter?

Minimal code demonstrating the bug or unexpected behavior

Minimal runnable code

using ITensors

N = 32
sites = siteinds("S=1/2", N; conserve_qns=false)

ampo = AutoMPO()
for j in 1:(N - 1)
  ampo .+= 0.5, "S+", j, "S-", j + 1
  ampo .+= 0.5, "S-", j, "S+", j + 1
  ampo .+= "Sz", j, "Sz", j + 1
end
H = MPO(ampo, sites)

ψ₀ = randomMPS(sites)

sweeps = Sweeps(6)
maxdim!(sweeps, 4, 16, 64)
mindim!(sweeps, 1, 16,64,)
cutoff!(sweeps,1e-10)
noise!(sweeps, 0)
energy, ψ = @time dmrg(H, ψ₀, sweeps)
@show linkdims(ψ)

Expected output or behavior

Expected Output

If cutoff!(sweeps,0):

linkdims(ψ) = [2, 4, 8, 16, 32, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 32, 16, 8, 4, 2]

Actual output or behavior

Output of minimal runnable code

linkdims(ψ) = [64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 32, 16, 8, 4, 2]

Version information

• Output from versioninfo():

julia> versioninfo()
Julia Version 1.8.5
Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 24 × Intel(R) Xeon(R) Gold 6128 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake-avx512)
  Threads: 1 on 24 virtual cores

• Output from using Pkg; Pkg.status("ITensors"):

julia> using Pkg; Pkg.status("ITensors")
Status `~/.julia/environments/v1.8/Project.toml`
  [9136182c] ITensors v0.3.43

[mtfishman]

I think this may be because we don't really have a full QR/SVD, only thin versions, so the maximum dimension is set by the smaller matrix dimension and you get these boundary effects depending on the sweep direction.

[mtfishman]

Though I'm a bit confused about what behavior you are hoping for, do you want the bonds to be 64 everywhere?

Because of the thin QR/SVD, the bond dimension pattern will depend on the gauge center, which isn't really something particular to the DMRG function (i.e. try orthogonalize(ψ, 1), orthogonalize(ψ, N), orthogonalize(ψ, N÷2)).

[ryanlevy]

My personal preference would be forcing of the bulk tensors to be the mindim but would also be ok with 64 everywhere, as long as its consistent. (This came up benchmarking two codes where I was trying to get the tensor sizes to match exactly)

[mtfishman]

I thought in your example above mindim is 64, so I'm a bit confused by your comment where it seems like you're implying they are different from each other. Also it seems like the bulk tensors are mindim/64 in the output you show in your first post.

[mtfishman]

To get the output you expected (based on your first post) does it work to call orthogonalize(ψ, N÷2)?

[ryanlevy]

Yeah I can call orthogonalize(ψ, N÷2) to get it in the right form at the end, but during DMRG sweeping the wave function sometimes has [2, 4, 8, 16, 32, 64, 64,... as the link dims but ends with [64, 64, 64, 64, 64, 64, 64, 64, on the left side only.
Here are the link dims at every sweep for N=12

[2, 64, 64, 64, 64, 64, 32, 16, 8, 4, 2]
[2, 4, 64, 64, 64, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 64, 64, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 16, 64, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 64, 16, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 64, 64, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 64, 64, 64, 4, 2]
[2, 4, 8, 16, 32, 64, 64, 64, 64, 64, 2]
[2, 4, 8, 16, 32, 64, 64, 64, 64, 64, 64]
[2, 4, 8, 16, 32, 64, 64, 64, 64, 64, 2]
[2, 4, 8, 16, 32, 64, 64, 64, 64, 4, 2]
[2, 4, 8, 16, 32, 64, 64, 64, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 64, 16, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 16, 32, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 16, 64, 64, 32, 16, 8, 4, 2]
[2, 4, 8, 64, 64, 64, 32, 16, 8, 4, 2]
[2, 4, 64, 64, 64, 64, 32, 16, 8, 4, 2]
[2, 64, 64, 64, 64, 64, 32, 16, 8, 4, 2]
[64, 64, 64, 64, 64, 64, 32, 16, 8, 4, 2]

[mtfishman]

Is this the bond dimension after each sweep or each 2-site update? If it's the latter, that's what I would expect based on how MPS gauging works.

[ryanlevy]

It's after each 2-site update (in measure! I have println("$(linkdims(psi))")).
After each sweep (once we hit 64) the bond dimensions are [64, 64, 64, 64, 64, 64, 32, 16, 8, 4, 2]

[mtfishman]

So I guess I'm not seeing an issue here, though happy to discuss offline why this is expected behavior.

[mtfishman]

Basically, if you imagine doing a sweep without doing any actual updates (say DMRG on an identity MPO), it should act like calling orthogonalize on each sweep location j.

[ryanlevy]

Quick update, a more minimal example. Say we have a two-site tensor that we want to factorize

s  = Index(2)
s′ = Index(2)
l = Index(64)
r = Index(16)
T = randomITensor(l,s,s′,r)

The inconsistent behavior manifests as

> L,R,spec = factorize(T,l,s; cutoff=0,mindim=64)
> @show inds(L)
inds(L) = ((dim=64|id=719), (dim=2|id=499), (dim=32|id=847|"Link,fact"))

where we get the 64->32->16 stepdown.
Or with cutoff

> L,R,spec = factorize(T,l,s; cutoff=1e-10,mindim=64)
> @show inds(L)
inds(L) = ((dim=64|id=444), (dim=2|id=941), (dim=64|id=961|"Link,fact"))

Where the factorize respects the mindim by padding the link dimension. I think the cutoff=0 behavior is best here.

[mtfishman]

Thanks. So when you call factorize with cutoff=1e-10, it calls the ITensors.factorize_eigen backend, which uses a Hermitian eigendecomposition instead of an SVD for performance. That backend is used automatically when the cutoff is larger that 1e-12 and therefore high precision of the singular values isn't needed.

You can see that here:

julia> L, R = ITensors.factorize_eigen(T, l, s; cutoff=1e-4, mindim=64);

julia> R
ITensor ord=3 (dim=64|id=630|"Link,eigen")' (dim=2|id=791) (dim=16|id=387)
NDTensors.Dense{Float64, Vector{Float64}}

julia> L
ITensor ord=3 (dim=64|id=131) (dim=2|id=103) (dim=64|id=630|"Link,eigen")'
NDTensors.Dense{Float64, Vector{Float64}}

This helps explain some of the confusion I had about this behavior, since using an eigendecomposition circumvents the thin QR/SVD behavior. It also explains why it occurs only for larger cutoff values.

I quick fix could be to set the maxdim in factorize_eigen to the minimum dimension of the left or right spaces.