Implement approx_itensornetwork

[LinjianMa]

Interfaces:

"""
Approximate a `binary_tree_partition` into an output ITensorNetwork
with the same binary tree structure. `root` is the root vertex of the
pre-order depth-first-search traversal used to perform the truncations.
"""
function approx_itensornetwork(
  ::Algorithm"density_matrix",
  binary_tree_partition::DataGraph;
  root,
  cutoff=1e-15,
  maxdim=10000,
  contraction_sequence_alg,
  contraction_sequence_kwargs,
)
"""
Approximate a given ITensorNetwork `tn` into an output ITensorNetwork with `output_structure`.
`output_structure` outputs a directed binary tree DataGraph defining the desired graph structure.
"""
function approx_itensornetwork(
  ::Algorithm"density_matrix",
  tn::ITensorNetwork,
  output_structure::Function=path_graph_structure;
  cutoff,
  maxdim,
  contraction_sequence_alg,
  contraction_sequence_kwargs,
)
"""
Approximate a given ITensorNetwork `tn` into an output ITensorNetwork
with a binary tree structure. The binary tree structure is defined based
on `inds_btree`, which is a directed binary tree DataGraph of indices.
"""
function approx_itensornetwork(
  ::Algorithm"density_matrix",
  tn::ITensorNetwork,
  inds_btree::DataGraph;
  cutoff=1e-15,
  maxdim=10000,
  contraction_sequence_alg="optimal",
  contraction_sequence_kwargs=(;),
)

Note: The current code is designed so that its performance is perfect (as I’m aware of, this should already be the state-of-the-art density matrix algorithm implementation, in terms of the caching implementations)

Some internally-used definitions:

density matrix:

Example

"""
Example:
  Consider a tensor network below,
    1
    /\
   9  2
  /   /\
     3  6
    /|  /\
   4 5 7  8
  /  | |   \
  The density matrix for the edge `NamedEdge(2, 3)` squares the subgraph with vertices 3, 4, 5
     |
     3
    /|
   4 5
   | |
   4 5
   |/
   3
   |
  The density matrix for the edge `NamedEdge(5, 3)` squares the subgraph
    with vertices 1, 2, 3, 4, 6, 7, 8, 9
      1
      /\
     /  2
    /   /\
   /   3  6
  9   /|  /\
  |  4   7  8
  |  |   |  |
  |  4   7  8
  |  |/  | /
  |  3   6
  |  |  /
  |  | /
  |  2
  9 /
  |/
  1
  The density matrix for the edge `NamedEdge(4, 3)` squares the subgraph
    with vertices 1, 2, 3, 5, 6, 7, 8, 9
      1
      /\
     /  2
    /   /\
   /   3  6
  9   /|  /\
  |    5 7  8
  |    | |  |
  |    5 7  8
  |  |/  | /
  |  3   6
  |  |  /
  |  | /
  |  2
  9 /
  |/
  1
"""

partial density matrix:

Example

"""
Example:
  Consider a tensor network below,
    1
    /\
   9  2
  /   /\
     3  6
    /|  /\
   4 5 7  8
  /  | |   \
  The partial density matrix for the Edge set `Set([NamedEdge(2, 3), NamedEdge(5, 3)])`
    squares the subgraph with vertices 4, and contract with the tensor 3
    |
    3
   /
  4 - 4 -
  The partial density matrix for the Edge set `Set([NamedEdge(4, 3), NamedEdge(5, 3)])`
    squares the subgraph with vertices 1, 2, 6, 7, 8, 9, and contract with the tensor 3
      1
      /\
     /  2
    /   /\
   /   3  6
  9   /|  /\
  |      7  8
  |      |  |
  |      7  8
  |      | /
  |      6
  |     /
  |  | /
  |  2
  9 /
  |/
  1
  The density matrix for the Edge set `Set([NamedEdge(4, 3), NamedEdge(2, 3)])`
    squares the subgraph with vertices 5. and contract with the tensor 3
    |
    3
   /
  5 - 5 -
"""

Calculate _DensityMatrix using _PartialDensityMatrix:

_PartialDensityMatrix is introduced for the efficiency purpose. For a density matrix dm with root v and children c1,c2, we can calculate dm via two _PartialDensityMatrix pdm1, pdm2, where pdm1 is defined on v with child c1 and pdm2 defined on v with child c2. The contraction can be done by property sim the inds of pdm2, then contract with pdm1. In this way, everytime when we want to use a specific dm we don’t need to re-compute it, since pdm1 or pdm2 can be cached and we can reuse them.

The density matrix algorithm core part is below:

 alg_graph = _DensityMartrixAlgGraph(partition, out_tree, root)
  output_tn = ITensorNetwork()
  for v in post_order_dfs_vertices(out_tree, root)[1:(end - 1)]
    U = _rem_vertex!(
      alg_graph, v; cutoff, maxdim, contraction_sequence_alg, contraction_sequence_kwargs
    )
    add_vertex!(output_tn, v)
    output_tn[v] = U
  end

each _rem_vertex! will remove v in the partition, and output its projector U. This U is then used to update the output tree tensor network. _rem_vertex! also update and use all the caches.

Discussions:
Note that it’s very possible that the _DensityMatrixAlgCaches has similar concepts as the ITensorNetworkCache @mtfishman plans to work on. If so I would be happy to discuss and we can think about how to design it in a better way.

[codecov-commenter]

Codecov Report

Merging #66 (fc5b85a) into main (536faaa) will increase coverage by 0.96%.
The diff coverage is 97.59%.

📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff             @@
##             main      #66      +/-   ##
==========================================
+ Coverage   76.82%   77.79%   +0.96%     
==========================================
  Files          60       61       +1     
  Lines        3115     3233     +118     
==========================================
+ Hits         2393     2515     +122     
+ Misses        722      718       -4     

Impacted Files	Coverage Δ
src/ITensorNetworks.jl	77.77% <ø> (ø)
src/contract_deltas.jl	86.76% <ø> (-0.74%)	⬇️
src/itensors.jl	82.85% <83.33%> (+0.09%)	⬆️
src/approx_itensornetwork.jl	97.91% <97.91%> (ø)
src/binary_tree_partition.jl	100.00% <100.00%> (ø)
src/partition.jl	86.13% <100.00%> (+1.52%)	⬆️

... and 9 files with indirect coverage changes

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[LinjianMa]

@mtfishman Not sure what happened with the documentation CI, any ideas?

[mtfishman]

Not sure, I hadn't even realized we had doctests. I don't think that doctest even makes sense since it relies on https://github.com/GiggleLiu/ITensorContractionOrders.jl which was moved into this package anyway, so we could just remove that doctest:
https://github.com/mtfishman/ITensorNetworks.jl/blob/64f17f1d7a567d2d6b7c56e551f84f2ff34a1f25/src/requires/omeinsumcontractionorders_itensorcontractiontree.jl#L12-L20

EDIT: It looks like that test starting getting run when it wasn't running before because you added using OMEinsumContractionOrders and added it as a dependency, which triggered that code to load when the package is loaded. I don't think we should do that, since it should be conditionally loaded externally by users (using https://github.com/JuliaPackaging/Requires.jl). But either way we should remove that doctest since it doesn't make sense to have in this package.

[mtfishman]

This makes me think about a slightly different interface and code structure. Here is a summary:

function approx_itensornetwork(tn::ITensorNetwork, output_structure::IndsNetwork; alg="density_matrix", cutoff, maxdim)
  partitioned_tn = partition(tn, output_structure) # Partition the network based on the `output_structure`, outputs a `DataGraph`
  return approx_itensornetwork(partitioned_tn; alg, cutoff, maxdim)
end

function approx_itensornetwork(partitioned_tn::DataGraph; alg="density_matrix", cutoff, maxdim)
  return approx_itensornetwork(Algorithm(alg), partitioned_tn; cutoff, maxdim)
end

function approx_itensornetwork(alg::Algorithm"density_matrix", partitioned_tn::DataGraph; cutoff, maxdim)
  @assert is_tree(partitioned_tn) # For now restrict the desired tensor network structure to be a tree
  # Implementation of the density matrix algorithm on the partitioned tensor network
end

function approx_itensornetwork(alg::Algorithm"orthogonalize", partitioned_tn::DataGraph; cutoff, maxdim)
  @assert is_tree(partitioned_tn) # For now restrict the desired tensor network structure to be a tree
  # Implementation of the orthogonalization algorithm on the partitioned tensor network.
  # Probably has overlap with the density matrix algorithm so code should be shared through generic functions.
end

function approx_itensornetwork(tn::ITensorNetwork, output_structure::Function=path_graph_structure; alg, cutoff, maxdim)
  output_structure_indsnetwork = output_structure(tn) # Outputs an `IndsNetwork`
  return approx_itensornetwork(tn, output_structure_indsnetwork; alg, cutoff, maxdim)
end

function path_graph_structure(tn::ITensorNetwork)
  # Outputs a maximimally unbalanced binary tree IndsNetwork defining the desired graph structure
end

function binary_tree_structure(tn::ITensorNetwork)
  # Outputs a binary tree IndsNetwork defining the desired graph structure
end

Besides some restructuring/renaming, I think the main difference would be using an IndsNetwork to define the desired structure of the output tensor network instead of a nested tree structure like you did in #42 and #64. The reasons would be:

1. I think an actual graph structure like IndsNetwork is easier to work with and easier for users to construct and understand.
2. It would anticipate generalizing the function to the case where the desired output network isn't a tree, but some more general tensor network (such as approximating the contraction of a 2D-like network as a PEPS/2D tensor network).

Let me know if that makes sense. Some of the proposal may clash with realities about the flow/logic of the algorithm so that would be helpful to know as well.

[LinjianMa]

@mtfishman The suggestions above as to the interfaces make sense to me, here are some details I would like to bring up:

1. In addition to the tree structure defined by indsnetwork, the density matrix algorithm and also the orthogonalize algorithm needs another input root to indicate the root vertex for the bfs traversal. We can add this in the function argument, but I'm not sure if this is always useful (for non-tree approximation we may not need this).
2. The current density matrix algorithm assumes that the bfs tree of indsnetwork with specific root is a binary tree. It's possible to generalize it to arbitrary trees, but it takes time for an efficient implementation (in particular, the caching logic needs to be designed carefully). So for now, I suggest we restrict the density matrix algorithm to only work with binary trees, and we could add a function like is_binary_tree to make sure the structure is allowed
3. As to partition(tn, output_structure): for now with the density matrix algorithm the partition function is called binary_tree_partition. We can make it a specialized partition function, for example

partition(alg::Algorithm"mincut_recursive_bisection", tn::ITensorNetwork, output_structure::IndsNetwork, root)

[mtfishman]

That's helpful, thanks. Maybe what we could do is use a directed DataGraph instead of an IndsNetwork, which would then have to be a binary arborescence/rooted directed tree in order for it to work with the current algorithm.

Definitely fine to keep it working only for binary trees, good idea to add a check for that. Also sounds like a good plan to specialize on partition.

[mtfishman]

Was there an issue using ITensorNetwork instead of ITensorNetwork{Any}, i.e. let it try to infer the vertex type?

[LinjianMa]

The reason I force it to be ITensorNetwork{Any} is that in approx_itensornetwork I update each partition inplace during the density matrix algorithm, and that update sometimes will assign a new itensornetwork type to each partition.

But it's right that forcing the type here doesn't make sense. I made a modification to update the datatype of partition inside approx_itensornetwork right before the update of density matrix algorithm, so the problem should be fixed now.

[LinjianMa]

@mtfishman could you retake a look at this?

[mtfishman]

@LinjianMa yes I was trying to make some time for this, I should have time in the next few days.

There is a lot of code in approx_itensornetwork.jl, can you write up a summary of the code logic in the first post, and include any subtle design choices that I should focus on?

Also, I have a more general question about how this fits into the over all algorithm. My current understanding is that the algorithm can be summarized as follows.

Algorithm
Contract the subgraphs Algorithm Input 2 of a tensor network Algorithm Input 1 according to the subgraph contraction sequence Algorithm Input 3, by traversing the contractions sequence and performing pairwise contractions of subgraphs using the Subalgorithm below.

Algorithm Inputs

• Algorithm Input 1: A tensor network we want to approximately contract.
• Algorithm Input 2: A partitioning of Algorithm Input 1 into subgraphs.
• Algorithm Input 3: A subgraph contraction sequence for the subgraphs Algorithm Input 2.

Subalgorithm
For each pairwise contraction of two subgraphs of Algorithm Input 2 according to the contraction sequence Algorithm Input 3, approximately contract them using the density matrix algorithm, outputting a binary tree that approximates the contraction of the two subgraphs.

Subalgorithm Inputs

• Subalgorithm Input 1: Two subgraphs from Algorithm Input 2, determined by traversing the contraction sequence Algorithm Input 3.
• Subalgorithm Input 2: A desired binary tree output structure (determined for example from the function binary_tree_structure introduced in Introduce binary_tree_structure, path_graph_structure, and rename binary_tree_partition #69, which takes as an input the two subgraphs Subalgorithm Input 1).

Questions

1. Is the function approx_itensornetwork introduce in this PR an implementation of the Subalgorithm? If so, am I correct to assume that Algorithm will be implemented in a future PR, as a simple function that traverses the contraction sequence tree Algorithm Input 3 and performs pairwise contractions of subgraphs from Algorithm Input 2 using the Subalgorithm (approx_itensornetwork)?
2. From previous meetings, my understanding was that the desired binary tree output structure (Subalgorithm Input 2) relied on information of the overall graph connectivity, i.e. made use of the overall graph structure of whatever part of original tensor network (Algorithm Input 1) remains uncontracted in order to determine the binary tree which helps minimize the number of swaps needed for future subgraph contractions through the Subalgorithm. That doesn't seem to be the case in the current code design, could you clarify that for me?

[LinjianMa]

@mtfishman Sure, I will summarize approx_itensornetwork.jl

For your questions:

1. Yes, this is correct.
2. Note that we have an interface

function approx_itensornetwork(
  tn::ITensorNetwork,
  inds_btree::DataGraph;
  alg::String,
  cutoff=1e-15,
  maxdim=10000,
  contraction_sequence_alg="optimal",
  contraction_sequence_kwargs=(;),
)

to let approx_itensornetwork be used in Algorithm Input 1, we will let inds_btree be a directed tree of indices decided by some other functions, and inds_btree will be chosen based on the original tensor network.

So the interface below

function approx_itensornetwork(
  tn::ITensorNetwork,
  output_structure::Function;
  alg::String,
  cutoff=1e-15,
  maxdim=10000,
  contraction_sequence_alg="optimal",
  contraction_sequence_kwargs=(;),
)

will not be used there since output_structure only knows the structure of tn. But this can still be useful so we keep this interface.

[mtfishman]

Ok, so to summarize: approx_itensornetwork can accept any binary tree structure, which could be chosen based on the rest of the tensor network but that logic isn't in this PR. So that will be some extra logic that will be implemented in a future PR, and approx_itensornetwork and binary_tree_structure (introduced in #69) are written generally enough so that they can be used once that logic is implemented. Is that correct? If so, that helps me understand the scope and context of this PR.

The part I'm still confused about is that it seems like in order to implement the logic of determining the output structure from the rest of the graph it would affect the function binary_tree_structure introduced in #69, which from what I can tell doesn't implement that logic.

[LinjianMa]

@mtfishman yes that's correct.

Note that the current binary_tree_structure and path_graph_structure functions are also general enough. The general API is

function binary_tree_structure(tn::ITensorNetwork, outinds::Vector{<:Index})

where outinds is part of noncommoninds of tn, and this can be used to get the tree ordering of part of the noncommoninds. This function and also some other logic can help us decide the swap-efficient binary tree structures.

[LinjianMa]

@mtfishman the first post is updated

[mtfishman]

Thanks for the summary @LinjianMa, that's helpful. I'll review this PR by the end of the week.

[LinjianMa]

@mtfishman The PR has been updated, please let me know if you have any other questions!

[mtfishman]

Looks great, thanks for all of the documentation and thanks for iterating on the code! This will be really useful.
I'll go ahead and merge.