[ITensors] Call out to TensorOperations.jl for optimal sequence finding (#1637)

Author: JoeyT1994

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensors"
 uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 authors = ["Matthew Fishman <[email protected]>", "Miles Stoudenmire <[email protected]>"]
-version = "0.8.10"
+version = "0.9.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -29,11 +29,13 @@ Zeros = "bd1ec220-6eb4-527a-9b49-e79c3db6233b"
 
 [weakdeps]
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [extensions]
 ITensorsHDF5Ext = "HDF5"
+ITensorsTensorOperationsExt = "TensorOperations"
 ITensorsVectorInterfaceExt = "VectorInterface"
 ITensorsZygoteRulesExt = "ZygoteRules"
 
@@ -58,6 +60,7 @@ SimpleTraits = "0.9.4"
 SparseArrays = "<0.0.1, 1.10"
 StaticArrays = "0.12, 1.0"
 Strided = "1.1, 2"
+TensorOperations = "5.1.4"
 TimerOutputs = "0.5.5"
 TupleTools = "1.2"
 VectorInterface = "0.4, 0.5"
@@ -68,5 +71,6 @@ julia = "1.10"
 [extras]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"

README.md

@@ -37,6 +37,8 @@ ITensors.jl is supported by the Flatiron Institute, a division of the Simons Fou
 
 ## News
 
+- March 26, 2025: ITensors.jl v0.9 has been released. This is a minor breaking change since the `optimal_contraction_sequence` function now passes to the `optimaltree` function from [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl). The `TensorOperations` package therefore needs to be loaded in order for `optimal_contraction_sequence` to be used or if the flag `ITensors.enable_contraction_sequence_optimization()` is switched on.
+
 - March 22, 2025: As part of the latest release of ITensors.jl (v0.8.3), all documentation related to MPS/MPO functionality has been moved to the [ITensorMPS.jl documentation](https://docs.itensor.org/ITensorMPS).
 
 - February 22, 2025: Please note that there were issues installing the latest version of ITensors.jl (ITensors.jl v0.8) in older versions of Julia v1.10 and v1.11 ([https://github.com/ITensor/ITensors.jl/issues/1618](https://github.com/ITensor/ITensors.jl/issues/1618), [https://itensor.discourse.group/t/typeparameteraccessors-not-found-error-on-julia-v-1-10-0/2260](https://itensor.discourse.group/t/typeparameteraccessors-not-found-error-on-julia-v-1-10-0/2260)). This issue has been fixed in [NDTensors.jl v0.4.4](https://github.com/ITensor/ITensors.jl/pull/1623), so please try updating your packages if you are using older versions of Julia v1.10 or v1.11 and running into issues installing ITensors.jl.

docs/Project.toml

@@ -4,10 +4,11 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
 Documenter = "1.9.0"
 HDF5 = "0.17.2"
-ITensors = "0.8.4"
 Strided = "2.2.0"
+TensorOperations = "5.1.4"
 LinearAlgebra = "1.10.0"

docs/src/ContractionSequenceOptimization.md

@@ -1,8 +1,6 @@
 # Contraction sequence optimization
 
-When contracting a tensor network, the sequence of contraction makes a big difference in the computational cost. However, the complexity of determining the optimal sequence grows exponentially with the number of tensors, but there are many heuristic algorithms available for computing optimal sequences for small networks[^1][^2][^3][^4][^5][^6]. ITensors.jl provides some functionality for helping you find the optimal contraction sequence for small tensor network, as we will show below.
-
-The algorithm in ITensors.jl currently uses a modified version of[^1] with simplifications for outer product contractions and is based on the implementation in [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl).
+When contracting a tensor network, the sequence of contraction makes a big difference in the computational cost. However, the complexity of determining the optimal sequence grows exponentially with the number of tensors, but there are many heuristic algorithms available for computing optimal sequences for small networks[^1][^2][^3][^4][^5][^6]. ITensors.jl imports functionality from [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl) for helping you find the optimal contraction sequence for small tensor network, as we will show below.
 
 [^1]: [Faster identification of optimal contraction sequences for tensor networks](https://arxiv.org/abs/1304.6112)
 [^2]: [Improving the efficiency of variational tensor network algorithms](https://arxiv.org/abs/1310.8023)
@@ -64,12 +62,13 @@ display(cost2)
 ```
 This example helps us learn that in the limit of large MPS bond dimension `m`, the first contraction sequence is faster, while in the limit of large MPO bond dimension `k`, the second sequence is faster. This has practical implications for writing an efficient DMRG algorithm in both limits, which we plan to incorporate into ITensors.jl.
 
-Here is a more systematic example of searching through the parameter space to find optimal contraction sequences:
+Here is a more systematic example of searching through the parameter space to find optimal contraction sequences. Note, the TensorOperations.jl library must be loaded to use the optimal_contraction_sequence function:
 ```julia
 using ITensors
 using Symbolics
 
 using ITensors: contraction_cost, optimal_contraction_sequence
+using TensorOperations: TensorOperations
 
 function tensor_network(; m, k, d)
   l = Index(m, "l")

ext/ITensorsTensorOperationsExt/ITensorsTensorOperationsExt.jl

@@ -0,0 +1,15 @@
+module ITensorsTensorOperationsExt
+
+using ITensors: ITensors, ITensor, dim, inds
+using TensorOperations: optimaltree
+
+function ITensors.optimal_contraction_sequence(
+  As::Union{Vector{<:ITensor},Tuple{Vararg{ITensor}}}
+)
+  network = collect.(inds.(As))
+  inds_to_dims = Dict(i => Float64(dim(i)) for i in unique(reduce(vcat, network)))
+  seq, _ = optimaltree(network, inds_to_dims)
+  return seq
+end
+
+end

src/ITensors.jl

@@ -51,9 +51,6 @@ Documentation: https://itensor.github.io/ITensors.jl/stable/
 module ITensors
 include("usings.jl")
 include("utils.jl")
-include("lib/ContractionSequenceOptimization/src/ContractionSequenceOptimization.jl")
-# TODO: `using .ContractionSequenceOptimization: ContractionSequenceOptimization, ...`.
-using .ContractionSequenceOptimization
 include("lib/LazyApply/src/LazyApply.jl")
 # TODO: `using .LazyApply: LazyApply, ...`.
 using .LazyApply
@@ -105,6 +102,7 @@ include("indexset.jl")
 include("itensor.jl")
 include("qn/flux.jl")
 include("oneitensor.jl")
+include("tensor_operations/contraction_cost.jl")
 include("tensor_operations/tensor_algebra.jl")
 include("tensor_operations/matrix_algebra.jl")
 include("tensor_operations/permutations.jl")

src/imports.jl

@@ -81,9 +81,6 @@ import Base.Broadcast:
   broadcastable,
   instantiate
 
-import ITensors.ContractionSequenceOptimization:
-  contraction_cost, optimal_contraction_sequence
-
 import Adapt: adapt_structure, adapt_storage
 
 import LinearAlgebra: