Current ortho fix

[JoeyT1994]

This PR fixes an issue with the current_ortho function which causes the alternating_update interface to break on more generic trees (current tests appeared to only check more standard cases like the comb tree so missed the bug). It also changes a test in test_dmrg to use a tree for which the code was previously breaking.

The issue appeared to be in the overlapping_vertex clause when at the end of a sweep which was assuming that the last vertex in the sweep was a neighbor of the previous vertex --- which may not always be true for the 'default_region_plan'. Removing the clause fixes the issue as the else clause always works anyway.

@mtfishman I feel like the "current_ortho" should either i) be a trait carried around by the state (by having, e.g., an orthogonal_tensornetwork type that wraps an ITensorNetwork and a v_orthocentre) or ii) something that the alternating update interface keeps track of as a variable and just re-assigns it each time a la

  psi = orthogonalize(psi, v, v_orthocentre)
  v_orthocentre = v

where orthogonalize can utilize the idea that if psi is a tree it just finds the shortest path between v and v_centre and does a sequence of QRs along that path for efficiency. If psi is not a tree it could run BP to put it in the Vidal gauge and then absorb square root lambdas everywhere except at v where it absorbs the neighboring lambda (making the BP environment identity at that site).

[JoeyT1994]

@mtfishman once you merge #207 I will update the formatting to be in line with main. I won't until then just to keep this PR easier to read.

[mtfishman]

@JoeyT1994 regarding the formatting, see ITensor/ITensors.jl#1536 (comment).

[mtfishman]

@mtfishman I feel like the "current_ortho" should either i) be a trait carried around by the state (by having, e.g., an orthogonal_tensornetwork type that wraps an ITensorNetwork and a v_orthocentre) or ii) something that the alternating update interface keeps track of as a variable and just re-assigns it each time a la

The TTN type is designed like i): https://github.com/ITensor/ITensorNetworks.jl/blob/v0.11.21/src/treetensornetworks/treetensornetwork.jl#L11-L12, do you mean you also want that for non-tree networks?

[JoeyT1994]

Ah I see okay I will leave the formatting as is then and avoid v2.0 of JuliaFormatter

[JoeyT1994]

@mtfishman I feel like the "current_ortho" should either i) be a trait carried around by the state (by having, e.g., an orthogonal_tensornetwork type that wraps an ITensorNetwork and a v_orthocentre) or ii) something that the alternating update interface keeps track of as a variable and just re-assigns it each time a la

The TTN type is designed like i): https://github.com/ITensor/ITensorNetworks.jl/blob/v0.11.21/src/treetensornetworks/treetensornetwork.jl#L11-L12, do you mean you also want that for non-tree networks?

I see I didn't know the TTN type had that. I guess I was just thinking more generally, we could have an OrthogonalITensornetwork <: ITensorNetwork type instead which carries a v_centre around (and if it's a tree it means the environment for that region is exactly identity, and if it's loopy it means the environment on that site/region is, under the BP approximation, identity). That fit nicely with the VidalITensorNetwork and ITensorNetwork types we already have and we know how to move between using BP (/ QRs if it's a tree).

If the TTN type carries the centre-site information what is the purpose of the current_ortho function here? Can't it just be pulled from the state the alternating_update routine is working with instead of the sweep plan?

[mtfishman]

If the TTN type carries the centre-site information what is the purpose of the current_ortho function here? Can't it just be pulled from the state the alternating_update routine is working with instead of the sweep plan?

That's a good question, I don't remember why it was designed that way. It looks like current_ortho was introduced by @b-kloss in #143, @b-kloss do you remember the logic behind that design vs. just using the orthogonality center information stored in the TTN struct?

[mtfishman]

@JoeyT1994 if it seems easier to just use the orthogonality center stored in the TTN, feel free to try to refactor it in that way.

[b-kloss]

If the TTN type carries the centre-site information what is the purpose of the current_ortho function here? Can't it just be pulled from the state the alternating_update routine is working with instead of the sweep plan?

That's a good question, I don't remember why it was designed that way. It looks like current_ortho was introduced by @b-kloss in #143, @b-kloss do you remember the logic behind that design vs. just using the orthogonality center information stored in the TTN struct?

The function current_ortho always returns a single vertex, while the orthogonality center of a TTN is a (contiguous on the graph) set of vertices. They are not equivalent, and the TTN by itself does not provide enough information to determine what the return of current_ortho would be inside alternating_update.

For nsite=1, i.e. single-site alternating update, current_ortho will return the orthogonality center.
For nsite=2, i.e. two-site alternating update, current_ortho will return the site that overlaps with the previous 2-site region in the sweep plan (which is the one which won't overlap with next 2-site region). For example, the default_inserter moves the orthogonality center to the vertex that's not current_ortho in the 2-site region. The logic inside current_ortho is convoluted because it becomes more complicated than what I described on general trees.

This logic may not generalize to nsite>2, and the whole mechanism can certainly be improved.

@JoeyT1994 The bugfix to current_ortho looks good to me and should work for generic trees. Would be good if there was a systematic way to test for generic trees, but that seems hard.

[JoeyT1994]

Thanks for the explanation @b-kloss !

So in the code current_ortho always returns a single vertex ortho_centre which is not actually the current orthogonality centre but is what you want the orthogonality centre to be? Then that current_ortho is used in a call to orthogonalize(ttn, ortho_centre) via the default_extracter?

Why does current_ortho need to be a single vertex anyway? Could we not just keep track of the ortho_region as the region on which current_ortho was last called (and thus it can be obtained as ortho_region(ttn)) and then during the extractor stage of the region_update we call

ttn = orthogonalize(ttn, region)

to orthogonalize onto the new region?

I'll give a refactor a go with this plan in mind as it should be possible to simplify the code along these lines, unless I'm missing something.

[JoeyT1994]

@mtfishman @b-kloss I have done a refactor of the code now. The changes done are as follows:

1. I refactored orthogonalize so that the functionality is almost entirely written in abstractitensornetwork.jl and is centered around a base function with takes a vector of edges and performs a corresponding sequence of QR decompositions.
2. The orthogonalize functionality in abstracttreetensornetwork.jl is now just a function orthogonalize(ttn::AbstractTreeTensorNetwork, region::Vector) which works out the sequence of edges needed to to pass to orthogonalize to move the orthogonality region of a ttn from where it currently is to region.
3. region_update now no longer needs to work out any orthogonality information via ortho_centre, it can just be pulled from the ttn directly when needed and orthogonalize(ttn, region) can just be called in the extracter. The ortho_centre function is now no longer needed and the bug this PR aimed to fix is gone.

I think this simplifies the logic of things quite a bit. Note, I also had to reverse the direction of the edges in the reverse sweep of TDVP to make things align in the code which I did by just adding a keyword argument as I find it rather tricky to navigate the code for generating sweep_plans in its current form.

[mtfishman]

@JoeyT1994 this definitely looks like a nice simplification of the code, though I'm having trouble picturing how the old code and new code are working diagrammatically. Let's meet to discuss just so I make sure I understand what's going on and can be convinced that the new code and old code are doing the same thing (besides the cases that the new code fixes).

This could help with something I was trying in the tensor network impurity solver. I was trying to automate finding the tree structure but for tree structures I was generating I was hitting a bug in the solver code, hopefully it was this bug that you're fixing here! (EDIT: I confirmed that this fixes the bug I was hitting.)

[JoeyT1994]

Thanks @mtfishman . I incorporated some of your comments and yes let's meet to discuss how the logic works here.

One slight flag here (which was also present in the previous version) is that the edge_sequence to move the ortho_region is computed via an O(L) time complexity function with L the number of tensors in the network.
Hence this actually makes the alternating_update sweeping on trees have an O(L^{2}) time complexity term (the sequence is computed every region_update). Although this is likely much smaller than the O(L Poly(chi)) scaling it is not really optimal and we should be able to write a more sophisticated function to move the ortho_region that is O(1) if the current ortho_region and new ortho_region are not extensive and are not separated by an extensive path.

[mtfishman]

Looks good, thanks!