Subspace expansion does not seem to properly choose cutoff

[ebertok]

I'm not sure whether this is expected behavior,
but when choosing a large cutoff for the subspace expansion (e.g. 10^-3), and a large maxdim,
the current subspace expansion will still max out the bond dimension until maxdim is reached, even though cutoff is already reached. So, subspace expansion will go beyond the cutoff of 10^-3 and report an S[end,end] = 10^-11.
This seems to happen mostly when outer_iters is large (I chose 20).

The way I understand the routine, it should automatically truncate down to the largest possible bond dimension that still suffices for the given cutoff.

[LHerviou]

Hi,
Is it possible that the contributions to S smaller than 10^-3 are already present in psi before the subspace expansion?
The current implementation (even before my modifications) does not control and truncate this part. The maxdim / cutoff only act on the added nullspace.

@mtfishman I do not know exactly what is the best intended behavior here.

[mtfishman]

Yes that's correct, the current implementation only expands the bond dimension, it never truncates to a smaller bond dimension than you already have, so if you already have smaller singular values than your cutoff than those would remain after the subspace expansion. Perhaps what is happening in your case is that the bond dimension gets saturated early on when the singular values are not small, but then later in the calculation after more iterations/sweeps the singular values get smaller but with the current limitations of the implementation they don't get truncated away.

Could you show the output of S[end,end] for your entire calculation?

We could add a separate step which truncates the infinite MPS according to the singular values. It simply involves doing the SVD of the C bond matrices as C = U*S*V, truncating according to the singular values S, and then changing the gauge of the left and right gauged MPS by U and V. However, this puts the MPS into the gauge where the bond matrices S are diagonal, and I seem to recall trying to run VUMPS in that gauge (many years ago) and it caused VUMPS to become numerically unstable. My understanding was that it was hard for the eigensolver to solve for eigenvectors which were diagonal in that way, or that it caused issues with the step of solving for AL and AR from AC and C. Perhaps a way to mitigate that problem would be to first truncate, and then rotate by a random unitary gauge so that the new bond matrices are no longer diagonal

Over all, I am hoping to make all of this more automated, where you just set a maxdim and cutoff like in the ITensors.jl dmrg function and then it automatically does the proper subspace expansion/truncation, but I haven't had time to work on that yet. The subspace expansion could be automatically enabled when the smallest singular values reach the cutoff (I came across a TDVP paper that used that technique but I can't find it now).