From 4994a212751fd0ef0dd91dfe11bac8005e6ffbfc Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 14:07:19 -0400
Subject: [PATCH 01/11] Delete insert_linkinds(::AbstractITensorNetwork, ...)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Migrate the two in-tree callers off `insert_linkinds` and onto
`add_edge!`:

- `add(tn1, tn2)` iterates the union of edges and calls `add_edge!` on
  each (no-op for edges already present).
- `tebd` does the same once before the time-evolution loop, replacing
  the per-step `insert_linkinds` call.

Strengthen `add_edge!` to also thread link inds when the graph edge
already exists but no shared Index is present (an inconsistent state
that today's API still allows). The long-term invariant is
graph-edge ↔ shared-index one-to-one; until that's enforced
structurally, `add_edge!` fills in whichever piece is missing. Return
value follows the `Graphs.add_edge!` convention: `true` only when a
new graph edge was added.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 src/abstractitensornetwork.jl | 36 ++++++++++++++---------------------
 src/tebd.jl                   |  5 ++++-
 2 files changed, 18 insertions(+), 23 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index ed98ed45..5f7baa00 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -824,21 +824,6 @@ function linkdims(tn::AbstractITensorNetwork{V}) where {V}
     return ld
 end
 
-function insert_linkinds(
-        tn::AbstractITensorNetwork, edges = edges(tn); link_space = trivial_space(tn)
-    )
-    tn = copy(tn)
-    for e in edges
-        if isempty(linkinds(tn, e))
-            iₑ = Index(link_space, edge_tag(e))
-            X = onehot(iₑ => 1)
-            tn[src(e)] *= X
-            tn[dst(e)] *= dag(X)
-        end
-    end
-    return tn
-end
-
 # Factor `tn[src(edge)]` via `f` along the inds it doesn't share with
 # `tn[dst(edge)]`, leaving the new factor on `src(edge)` and contracting
 # the residual into `tn[dst(edge)]`. With `f = qr` and no shared inds yet,
@@ -853,11 +838,17 @@ function factorize_edge!(f, tn::AbstractITensorNetwork, edge)
     return tn
 end
 
+# Ensure `edge` is present in both the underlying graph and as a shared
+# Index between the endpoint tensors. The long-term invariant is that the
+# two are one-to-one; today they can drift apart, so this fills in
+# whichever piece is missing. Returns the `Graphs.add_edge!` convention:
+# `true` if a new graph edge was added, `false` if the graph edge was
+# already present (regardless of whether link inds had to be threaded).
 function Graphs.add_edge!(tn::AbstractITensorNetwork, edge)
-    has_edge(tn, edge) && return false
-    add_edge!(data_graph(tn), edge)
-    factorize_edge!(qr, tn, edge)
-    return true
+    added = !has_edge(tn, edge)
+    added && add_edge!(data_graph(tn), edge)
+    isempty(linkinds(tn, edge)) && factorize_edge!(qr, tn, edge)
+    return added
 end
 
 # TODO: What to output? Could be an `IndsNetwork`. Or maybe
@@ -900,9 +891,10 @@ function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
 
     if !issetequal(edges_tn1, edges_tn2)
-        new_edges = union(edges_tn1, edges_tn2)
-        tn1 = insert_linkinds(tn1, new_edges)
-        tn2 = insert_linkinds(tn2, new_edges)
+        for e in union(edges_tn1, edges_tn2)
+            add_edge!(tn1, e)
+            add_edge!(tn2, e)
+        end
     end
 
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
diff --git a/src/tebd.jl b/src/tebd.jl
index 5f4891f1..592fc852 100644
--- a/src/tebd.jl
+++ b/src/tebd.jl
@@ -15,11 +15,14 @@ function tebd(
     s = siteinds(ψ)
     u⃗ = Vector{ITensor}(𝒰, s)
     nsteps = Int(β ÷ Δβ)
+    ψ = copy(ψ)
+    for e in edges(ψ)
+        add_edge!(ψ, e)
+    end
     for step in 1:nsteps
         if step % print_frequency == 0
             @show step, (step - 1) * Δβ, β
         end
-        ψ = insert_linkinds(ψ)
         ψ = apply(u⃗, ψ; cutoff, maxdim, normalize = true, ortho, kwargs...)
         if ortho
             for v in vertices(ψ)

From 86adbfa7c0be5c91f9a81f6ca2c5c2219ef2aa18 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 14:10:30 -0400
Subject: [PATCH 02/11] Use GraphsExtensions.add_edges! for the multi-edge
 migration sites

Replaces the manual per-edge loops in `add(tn1, tn2)` and the tebd
initial-state setup with `add_edges!`, which dispatches per-edge to
the reformed `add_edge!` (ensures graph edge + link inds).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 src/abstractitensornetwork.jl | 9 ++++-----
 src/tebd.jl                   | 4 +---
 2 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 5f7baa00..5c6d7cca 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -12,7 +12,7 @@ using LinearAlgebra: LinearAlgebra, factorize, qr
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
-    directed_graph, incident_edges, rename_vertices, vertextype, ⊔
+    add_edges!, directed_graph, incident_edges, rename_vertices, vertextype, ⊔
 using NamedGraphs: NamedGraphs, NamedGraph, Vertices, not_implemented, steiner_tree
 using SplitApplyCombine: flatten
 
@@ -891,10 +891,9 @@ function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
 
     if !issetequal(edges_tn1, edges_tn2)
-        for e in union(edges_tn1, edges_tn2)
-            add_edge!(tn1, e)
-            add_edge!(tn2, e)
-        end
+        new_edges = union(edges_tn1, edges_tn2)
+        add_edges!(tn1, new_edges)
+        add_edges!(tn2, new_edges)
     end
 
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
diff --git a/src/tebd.jl b/src/tebd.jl
index 592fc852..56ac99af 100644
--- a/src/tebd.jl
+++ b/src/tebd.jl
@@ -16,9 +16,7 @@ function tebd(
     u⃗ = Vector{ITensor}(𝒰, s)
     nsteps = Int(β ÷ Δβ)
     ψ = copy(ψ)
-    for e in edges(ψ)
-        add_edge!(ψ, e)
-    end
+    add_edges!(ψ, edges(ψ))
     for step in 1:nsteps
         if step % print_frequency == 0
             @show step, (step - 1) * Δβ, β

From 414e2b62fd009104ec96242df67b3840bd83f79d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 14:24:59 -0400
Subject: [PATCH 03/11] Delete insert_linkinds(::AbstractIndsNetwork, ...)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Rewrite the only caller, `TreeTensorNetwork(::ITensor, ::IndsNetwork)`,
to build up the vertex-tensor `Dict` directly during the post-order DFS
instead of scaffolding a placeholder `ITensorNetwork` upfront. The
scaffold tensors' link inds (filled in by `insert_linkinds`) never
participated in the factorization — they were stripped by
`setdiff(inds(a), inds(tn[dst(e)]))` and then overwritten by the
`tn[src(e)] = a_l` writes anyway. The literal effective behavior was
`left_inds = setdiff(inds(a), site_inds[dst(e)])`, which is what the
new code computes directly.

With the IndsNetwork-form caller gone, the function (and its
`experimental_methods.md` entry) follow.

This does not change the algorithm's observable behavior — including
the existing quirk that the resulting TTN does not put each
`site_inds[v]` on its own vertex (the first vertex processed in
post-order absorbs every site ind that isn't `dst(e)`'s, etc.). That
structural issue is real but out of scope here; tracked separately.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 docs/src/experimental_methods.md            |  6 ------
 src/abstractindsnetwork.jl                  | 21 --------------------
 src/treetensornetworks/treetensornetwork.jl | 22 ++++++---------------
 3 files changed, 6 insertions(+), 43 deletions(-)

diff --git a/docs/src/experimental_methods.md b/docs/src/experimental_methods.md
index 0b7ca2ad..4170377e 100644
--- a/docs/src/experimental_methods.md
+++ b/docs/src/experimental_methods.md
@@ -224,12 +224,6 @@ Methods which still need to be discussed, modified, or deprecated.
   union_all_inds(is_in::AbstractIndsNetwork...)
   ```
 
-* Insert a default link index on every edge of `indsnetwork` that doesn't already
-  have one — `link_space` controls the default bond dimension (`abstractindsnetwork.jl`):
-  ```julia
-  insert_linkinds(indsnetwork::AbstractIndsNetwork, edges = edges(indsnetwork); link_space = trivial_space(indsnetwork))
-  ```
-
 #### IndsNetwork
 
 * Type-parameter accessors and graph-type metadata (`indsnetwork.jl`):
diff --git a/src/abstractindsnetwork.jl b/src/abstractindsnetwork.jl
index 49af846e..f00b2a56 100644
--- a/src/abstractindsnetwork.jl
+++ b/src/abstractindsnetwork.jl
@@ -93,24 +93,3 @@ function union_all_inds(is_in::AbstractIndsNetwork...)
     end
     return is_out
 end
-
-function insert_linkinds(
-        indsnetwork::AbstractIndsNetwork,
-        edges = edges(indsnetwork);
-        link_space = trivial_space(indsnetwork)
-    )
-    indsnetwork = copy(indsnetwork)
-    for e in edges
-        # TODO: Change to check if it is empty.
-        if !isassigned(indsnetwork, e)
-            if !isnothing(link_space)
-                iₑ = Index(link_space, edge_tag(e))
-                # TODO: Allow setting with just `Index`.
-                indsnetwork[e] = [iₑ]
-            else
-                indsnetwork[e] = []
-            end
-        end
-    end
-    return indsnetwork
-end
diff --git a/src/treetensornetworks/treetensornetwork.jl b/src/treetensornetworks/treetensornetwork.jl
index 7dc7cb25..8b487453 100644
--- a/src/treetensornetworks/treetensornetwork.jl
+++ b/src/treetensornetworks/treetensornetwork.jl
@@ -147,25 +147,15 @@ function TreeTensorNetwork(
         @assert hasinds(a, is[v])
     end
     @assert ortho_region ⊆ vertices(is)
-    is = insert_linkinds(is)
     g = NamedGraph(underlying_graph(is))
-    ts = Dict{vertextype(g), ITensor}()
-    for v in vertices(g)
-        site_inds = get(is, v, Index[])
-        edges_v = [edgetype(is)(v, nv) for nv in neighbors(is, v)]
-        link_inds = reduce(vcat, (is[e] for e in edges_v); init = Index[])
-        ts[v] = ITensor(site_inds..., link_inds...)
-    end
-    tn = ITensorNetwork(ts, g)
     ortho_center = first(ortho_region)
-    for e in post_order_dfs_edges(tn, ortho_center)
-        left_inds = setdiff(inds(a), inds(tn[dst(e)]))
+    ts = Dict{vertextype(g), ITensor}()
+    for e in post_order_dfs_edges(g, ortho_center)
+        left_inds = setdiff(inds(a), get(is, dst(e), Index[]))
         a_l, a_r = factorize(a, left_inds; tags = edge_tag(e), ortho = "left", kwargs...)
-        tn[src(e)] = a_l
-        is[e] = commoninds(a_l, a_r)
+        ts[src(e)] = a_l
         a = a_r
     end
-    tn[ortho_center] = a
-    ttn_a = TreeTensorNetwork(tn)
-    return orthogonalize(ttn_a, ortho_center)
+    ts[ortho_center] = a
+    return orthogonalize(TreeTensorNetwork(ITensorNetwork(ts, g)), ortho_center)
 end

From f704f6607c107f0b388590e741519a7f42c67b8b Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 14:51:36 -0400
Subject: [PATCH 04/11] Delete TreeTensorNetwork(::ITensor, ::IndsNetwork)
 constructor
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The constructor had no in-src callers and only two test-only callers
(`test/test_ttns.jl`'s 'Construct TTN from ITensor or Array' and
`test/test_ttno.jl`'s 'Construct TTN operator from ITensor or Array'),
each of which only verified that `contract(ttn) ≈ original_dense`.
Tracing the algorithm showed it was producing factorizations that
contract back to the original but do not lay out site indices on their
designated vertices — the first vertex processed in post-order absorbs
every site index that isn't `dst(e)`'s, leaving later vertices without
their own site indices. The round-trip tests passed despite this
because the misplaced site indices and trivial links still contract
correctly.

Drop the constructor and its two test callers rather than rewrite the
algorithm: nothing outside the tests was relying on it, and a correct
replacement can be designed if a real use case appears.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 src/treetensornetworks/treetensornetwork.jl | 41 +--------------------
 test/test_ttno.jl                           | 38 -------------------
 test/test_ttns.jl                           | 29 +--------------
 3 files changed, 3 insertions(+), 105 deletions(-)
 delete mode 100644 test/test_ttno.jl

diff --git a/src/treetensornetworks/treetensornetwork.jl b/src/treetensornetworks/treetensornetwork.jl
index 8b487453..01588063 100644
--- a/src/treetensornetworks/treetensornetwork.jl
+++ b/src/treetensornetworks/treetensornetwork.jl
@@ -1,10 +1,6 @@
-using DataGraphs: underlying_graph
 using Dictionaries: Indices
-using Graphs: path_graph
-using ITensors: ITensor
-using LinearAlgebra: factorize, normalize
-using NamedGraphs.GraphsExtensions: GraphsExtensions, vertextype
-using NamedGraphs: NamedGraph, similar_graph
+using NamedGraphs.GraphsExtensions: vertextype
+using NamedGraphs: similar_graph
 
 """
     TreeTensorNetwork{V} <: AbstractTreeTensorNetwork{V}
@@ -126,36 +122,3 @@ end
 function set_ortho_region(tn::TTN, ortho_region)
     return TreeTensorNetwork(tn.tensornetwork; ortho_region)
 end
-
-"""
-    TreeTensorNetwork(a::ITensor, is::IndsNetwork; ortho_region=..., kwargs...) -> TreeTensorNetwork
-
-Decompose a dense `ITensor` `a` into a `TreeTensorNetwork` with the tree structure
-described by the `IndsNetwork` `is`.
-
-Successive QR/SVD factorizations are applied following a post-order DFS traversal from the
-root vertex, then the network is orthogonalized to `ortho_region` (defaults to the root).
-Extra `kwargs` (e.g. `cutoff`, `maxdim`) are forwarded to the factorization.
-"""
-function TreeTensorNetwork(
-        a::ITensor,
-        is::IndsNetwork;
-        ortho_region = Indices([GraphsExtensions.default_root_vertex(is)]),
-        kwargs...
-    )
-    for v in vertices(is)
-        @assert hasinds(a, is[v])
-    end
-    @assert ortho_region ⊆ vertices(is)
-    g = NamedGraph(underlying_graph(is))
-    ortho_center = first(ortho_region)
-    ts = Dict{vertextype(g), ITensor}()
-    for e in post_order_dfs_edges(g, ortho_center)
-        left_inds = setdiff(inds(a), get(is, dst(e), Index[]))
-        a_l, a_r = factorize(a, left_inds; tags = edge_tag(e), ortho = "left", kwargs...)
-        ts[src(e)] = a_l
-        a = a_r
-    end
-    ts[ortho_center] = a
-    return orthogonalize(TreeTensorNetwork(ITensorNetwork(ts, g)), ortho_center)
-end
diff --git a/test/test_ttno.jl b/test/test_ttno.jl
deleted file mode 100644
index c876a3d9..00000000
--- a/test/test_ttno.jl
+++ /dev/null
@@ -1,38 +0,0 @@
-using Graphs: vertices
-using ITensorNetworks: TreeTensorNetwork, contract, ortho_region, siteinds, union_all_inds
-using ITensors: @disable_warn_order, prime, random_itensor
-using LinearAlgebra: norm
-using NamedGraphs.NamedGraphGenerators: named_comb_tree
-using Random: shuffle
-using StableRNGs: StableRNG
-using Test: @test, @testset
-@testset "TTN operator Basics" begin
-    # random comb tree
-    rng = StableRNG(1234)
-    tooth_lengths = rand(rng, 1:3, rand(rng, 2:4))
-    c = named_comb_tree(tooth_lengths)
-    # specify random site dimension on every site
-    dmap = v -> rand(rng, 1:3)
-    is = siteinds(dmap, c)
-    # operator site inds
-    is_isp = union_all_inds(is, prime(is; links = []))
-    # specify random linear vertex ordering of graph vertices
-    vertex_order = shuffle(rng, collect(vertices(c)))
-    @testset "Construct TTN operator from ITensor or Array" begin
-        cutoff = 1.0e-10
-        sites_o = [is_isp[v] for v in vertex_order]
-        # create random ITensor with these indices
-        rng = StableRNG(1234)
-        O = random_itensor(rng, sites_o...)
-        # dense TTN constructor from IndsNetwork
-        @disable_warn_order o1 = TreeTensorNetwork(O, is_isp; cutoff)
-        root_vertex = only(ortho_region(o1))
-        @disable_warn_order begin
-            O1 = contract(o1, root_vertex)
-        end
-        @test norm(O - O1) < 1.0e2 * cutoff
-    end
-    @testset "Ortho" begin
-        # TODO
-    end
-end
diff --git a/test/test_ttns.jl b/test/test_ttns.jl
index 64fac579..46278013 100644
--- a/test/test_ttns.jl
+++ b/test/test_ttns.jl
@@ -1,39 +1,12 @@
 using DataGraphs: vertex_data
 using Graphs: vertices
 using ITensorNetworks:
-    ITensorNetwork, TreeTensorNetwork, contract, ortho_region, orthogonalize, siteinds
-using ITensors: @disable_warn_order, random_itensor
-using LinearAlgebra: norm
+    ITensorNetwork, TreeTensorNetwork, ortho_region, orthogonalize, siteinds
 using NamedGraphs.NamedGraphGenerators: named_comb_tree
-using Random: shuffle
 using StableRNGs: StableRNG
 using Test: @test, @testset
 include("utils.jl")
 @testset "TTN Basics" begin
-    # random comb tree
-    rng = StableRNG(1234)
-    tooth_lengths = rand(rng, 1:3, rand(rng, 2:4))
-    c = named_comb_tree(tooth_lengths)
-    # specify random site dimension on every site
-    dmap = v -> rand(rng, 1:3)
-    is = siteinds(dmap, c)
-    # specify random linear vertex ordering of graph vertices
-    vertex_order = shuffle(rng, collect(vertices(c)))
-
-    @testset "Construct TTN from ITensor or Array" begin
-        cutoff = 1.0e-10
-        # create random ITensor with these indices
-        rng = StableRNG(1234)
-        S = random_itensor(rng, vertex_data(is)...)
-        # dense TTN constructor from IndsNetwork
-        @disable_warn_order s1 = TreeTensorNetwork(S, is; cutoff)
-        root_vertex = only(ortho_region(s1))
-        @disable_warn_order begin
-            S1 = contract(s1, root_vertex)
-        end
-        @test norm(S - S1) < 1.0e2 * cutoff
-    end
-
     @testset "Convert ITN <-> TTN" begin
         g = named_comb_tree((3, 2))
         sites = siteinds("S=1/2", g)

From 6f9c31f35910634cdd76659c772676d027d5a01d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 15:00:14 -0400
Subject: [PATCH 05/11] Drop docs/src/tree_tensor_networks.md section for the
 deleted ITensor ctor

The 'From a dense ITensor' section documented
`TreeTensorNetwork(::ITensor, ::IndsNetwork)`, which was deleted in
f704f66. The `@example` block would have crashed the doc build.
---
 docs/src/tree_tensor_networks.md | 19 -------------------
 1 file changed, 19 deletions(-)

diff --git a/docs/src/tree_tensor_networks.md b/docs/src/tree_tensor_networks.md
index ff740ffb..9f2411cf 100644
--- a/docs/src/tree_tensor_networks.md
+++ b/docs/src/tree_tensor_networks.md
@@ -68,25 +68,6 @@ ITensorNetworks.TreeTensorNetwork
 ITensorNetworks.ITensorNetwork(::ITensorNetworks.TreeTensorNetwork)
 ```
 
-### From a dense `ITensor`
-
-A dense tensor can be decomposed into a TTN by successive QR/SVD factorisations along the
-tree edges. Truncation parameters (e.g. `cutoff`, `maxdim`) are forwarded to the
-factorisation step.
-
-```@example main
-g_small = named_comb_tree((3, 1))
-sites_small = siteinds("S=1/2", g_small)
-A = ITensors.random_itensor(
-    only(sites_small[(1, 1)]), only(sites_small[(2, 1)]), only(sites_small[(3, 1)])
-)
-ttn_A = TreeTensorNetwork(A, sites_small)
-```
-
-```@docs; canonical=false
-ITensorNetworks.TreeTensorNetwork(::ITensors.ITensor, ::ITensorNetworks.IndsNetwork)
-```
-
 ## Orthogonal Gauge
 
 One of the most powerful features of tree tensor networks is the ability to bring the

From f1b74c3a10c948eac8827baed3b542c0efc9883e Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 15:19:07 -0400
Subject: [PATCH 06/11] Use the non-mutating GraphsExtensions.add_edges
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`add(tn1, tn2)` and `tebd` were both calling `add_edges!` on local
copies (already returned by `combine_linkinds` or a preceding
`copy(ψ)`). Use `add_edges` (non-mutating) directly: same semantics,
clearer functional style, and the explicit `copy(ψ)` in `tebd` is no
longer needed.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 src/abstractitensornetwork.jl | 6 +++---
 src/tebd.jl                   | 3 +--
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 5c6d7cca..1f6c417f 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -12,7 +12,7 @@ using LinearAlgebra: LinearAlgebra, factorize, qr
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
-    add_edges!, directed_graph, incident_edges, rename_vertices, vertextype, ⊔
+    add_edges, directed_graph, incident_edges, rename_vertices, vertextype, ⊔
 using NamedGraphs: NamedGraphs, NamedGraph, Vertices, not_implemented, steiner_tree
 using SplitApplyCombine: flatten
 
@@ -892,8 +892,8 @@ function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
 
     if !issetequal(edges_tn1, edges_tn2)
         new_edges = union(edges_tn1, edges_tn2)
-        add_edges!(tn1, new_edges)
-        add_edges!(tn2, new_edges)
+        tn1 = add_edges(tn1, new_edges)
+        tn2 = add_edges(tn2, new_edges)
     end
 
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
diff --git a/src/tebd.jl b/src/tebd.jl
index 56ac99af..a9de4d00 100644
--- a/src/tebd.jl
+++ b/src/tebd.jl
@@ -15,8 +15,7 @@ function tebd(
     s = siteinds(ψ)
     u⃗ = Vector{ITensor}(𝒰, s)
     nsteps = Int(β ÷ Δβ)
-    ψ = copy(ψ)
-    add_edges!(ψ, edges(ψ))
+    ψ = add_edges(ψ, edges(ψ))
     for step in 1:nsteps
         if step % print_frequency == 0
             @show step, (step - 1) * Δβ, β

From b9fa3517ff03e745c7dd2be0e8d9c03694e2c905 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 15:26:54 -0400
Subject: [PATCH 07/11] Delete combine_linkinds / linkinds_combiners /
 is_multi_edge
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The only caller — the multi-edge collapse step inside `add(tn1, tn2)` —
is replaced by an inline loop calling `factorize_edge!(qr, tn, e)` on
multi-edges. `factorize_edge!(qr, ...)` produces the same single-shared-
Index end state as `combine_linkinds`, plus a free rank truncation to
the actual rank of `tn[src]` on the (left-inds | multi-link)
bipartition. With no callers left, drop `combine_linkinds`,
`linkinds_combiners`, and `is_multi_edge`, and remove the
`combine_linkinds` entries from `docs/src/experimental_methods.md`.
The `combiner` import in `abstractitensornetwork.jl` is no longer
needed either.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 docs/src/experimental_methods.md |  9 ------
 src/abstractitensornetwork.jl    | 54 ++++++++------------------------
 2 files changed, 13 insertions(+), 50 deletions(-)

diff --git a/docs/src/experimental_methods.md b/docs/src/experimental_methods.md
index 4170377e..de602acc 100644
--- a/docs/src/experimental_methods.md
+++ b/docs/src/experimental_methods.md
@@ -4,15 +4,6 @@ Methods which still need to be discussed, modified, or deprecated.
 
 ## ITensorNetwork Methods
 
-* Combine (fuse) every link index of a tensor network, or a chosen set of edges, into
-  a single index per edge using `combiner` tensors. (`abstractitensornetwork.jl`).
-  Comment: it may be better to fold this into a more general interface `gauge_transform(f, tn::AbstractITensorNetwork)` that applies a gauge transformation $X\_e$ $X\_e^{-1}$ on each edge $e$.
-  ```julia
-  linkinds_combiners(tn::AbstractITensorNetwork; edges = edges(tn))
-  combine_linkinds(tn::AbstractITensorNetwork, combiners)
-  combine_linkinds(tn::AbstractITensorNetwork; edges = edges(tn))
-  ```
-
 * Extract the `IndsNetwork` of a tensor network — site indices per vertex and link
   indices per edge (`abstractitensornetwork.jl`):
   ```julia
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 1f6c417f..a0ac6814 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -5,9 +5,9 @@ using DataGraphs: DataGraphs, edge_data, get_vertex_data, is_vertex_assigned,
 using Dictionaries: Dictionary
 using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edges, edgetype,
     has_edge, ne, neighbors, rem_edge!, src, vertices
-using ITensors: ITensors, @Algorithm_str, ITensor, addtags, combiner, commoninds,
-    commontags, contract, dag, inds, noprime, onehot, prime, replaceprime, replacetags,
-    setprime, settags, sim, swaptags, tags
+using ITensors: ITensors, @Algorithm_str, ITensor, addtags, commoninds, commontags,
+    contract, dag, inds, noprime, onehot, prime, replaceprime, replacetags, setprime,
+    settags, sim, swaptags, tags
 using LinearAlgebra: LinearAlgebra, factorize, qr
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
@@ -709,38 +709,6 @@ function neighbor_vertices(ψ::AbstractITensorNetwork, T::ITensor)
     return first.(v⃗)
 end
 
-function linkinds_combiners(tn::AbstractITensorNetwork; edges = edges(tn))
-    combiners = DataGraph(
-        directed_graph(underlying_graph(tn));
-        vertex_data_type = ITensor,
-        edge_data_type = ITensor
-    )
-    for e in edges
-        C = combiner(linkinds(tn, e); tags = edge_tag(e))
-        combiners[e] = C
-        combiners[reverse(e)] = dag(C)
-    end
-    return combiners
-end
-
-function combine_linkinds(tn::AbstractITensorNetwork, combiners)
-    combined_tn = copy(tn)
-    for e in edges(tn)
-        if !isempty(linkinds(tn, e)) && isassigned(combiners, e)
-            combined_tn[src(e)] = combined_tn[src(e)] * combiners[e]
-            combined_tn[dst(e)] = combined_tn[dst(e)] * combiners[reverse(e)]
-        end
-    end
-    return combined_tn
-end
-
-function combine_linkinds(
-        tn::AbstractITensorNetwork; edges::Vector{<:Union{Pair, AbstractEdge}} = edges(tn)
-    )
-    combiners = linkinds_combiners(tn; edges)
-    return combine_linkinds(tn, combiners)
-end
-
 function inner_network(x::AbstractITensorNetwork, y::AbstractITensorNetwork; kwargs...)
     return LinearFormNetwork(x, y; kwargs...)
 end
@@ -865,10 +833,6 @@ function ITensors.commoninds(tn1::AbstractITensorNetwork, tn2::AbstractITensorNe
     return inds
 end
 
-# Check if the edge of an itensornetwork has multiple indices
-is_multi_edge(tn::AbstractITensorNetwork, e) = length(linkinds(tn, e)) > 1
-is_multi_edge(tn::AbstractITensorNetwork) = Base.Fix1(is_multi_edge, tn)
-
 """
     add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork) -> ITensorNetwork
 
@@ -885,8 +849,16 @@ See also: `Base.:+` for `TreeTensorNetwork`, `truncate`.
 function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
     @assert issetequal(vertices(tn1), vertices(tn2))
 
-    tn1 = combine_linkinds(tn1; edges = filter(is_multi_edge(tn1), edges(tn1)))
-    tn2 = combine_linkinds(tn2; edges = filter(is_multi_edge(tn2), edges(tn2)))
+    # Collapse any multi-edges (edges carrying >1 shared Index) to a single
+    # shared Index, so the per-edge direct sum below sees one link per edge.
+    tn1 = copy(tn1)
+    tn2 = copy(tn2)
+    for e in edges(tn1)
+        length(linkinds(tn1, e)) > 1 && factorize_edge!(qr, tn1, e)
+    end
+    for e in edges(tn2)
+        length(linkinds(tn2, e)) > 1 && factorize_edge!(qr, tn2, e)
+    end
 
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
 

From 6ce01760b2cdbec61bb59fcc4b8c4bd07064d98c Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 15:51:17 -0400
Subject: [PATCH 08/11] Consolidate edge factorization on qr / qr!

Replace the family of vertex-splitting edge factorizations
(`svd(tn, edge)`, `factorize(tn, edge)`, `qr(tn, edge)`) and the
internal `gauge_edge` + `factorize_edge!` helpers with a single
absorbing-into-dst `qr` / `qr!` pair:

- `qr!(tn, edge)` mutates `tn` in place: factorize `tn[src(edge)]`
  along the inds it doesn't share with `tn[dst(edge)]`, write the
  orthogonal factor `Q` to `src(edge)`, and absorb the residual `R`
  into `tn[dst(edge)]`. The vertex set is unchanged; the bond between
  the two endpoints is re-gauged.
- `qr(tn, edge) = qr!(copy(tn), edge)` is the non-mutating wrapper.

Drop the dispatch on `alg::Algorithm` from `gauge_walk`, `tree_gauge`,
and `gauge`. There was only ever one implementation
(`Algorithm"orthogonalize"`), so the parameter was vestigial. Public
`orthogonalize` / `tree_orthogonalize` signatures are unchanged.

Migrate callers:
- `add_edge!` and `add(tn1, tn2)` now call `qr!` directly instead of
  the deleted `factorize_edge!(qr, tn, edge)`.
- `gauge_walk` body now calls `qr!` directly instead of `gauge_edge`.

Update `docs/src/interface_methods.md` to drop the entries for
`svd(tn, edge)` and `factorize(tn, edge)` and to describe the new
absorbing semantics of `qr(tn, edge)`.

Delete the now-orphaned `factorize(tn, 4 => 3)` regression in
`test/test_itensornetwork.jl` (it was testing the vertex-splitting
behavior, asserting `nv(tn_ortho) == nv(tn) + 1`). The remaining
`tree_orthogonalize` checks in the same testset cover the absorbing
path.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 docs/src/interface_methods.md                 |  23 +--
 src/abstractitensornetwork.jl                 | 194 ++++--------------
 .../abstracttreetensornetwork.jl              |  18 +-
 test/test_itensornetwork.jl                   |   7 -
 4 files changed, 49 insertions(+), 193 deletions(-)

diff --git a/docs/src/interface_methods.md b/docs/src/interface_methods.md
index 7eb40df5..b628ea07 100644
--- a/docs/src/interface_methods.md
+++ b/docs/src/interface_methods.md
@@ -48,21 +48,14 @@ These ITensorNetwork constructor interfaces are foundational to other constructo
   contract(tn::AbstractITensorNetwork, edge::Pair; kws...)
   ```
 
-* Factorize the bond on `edge` using the default factorization (`abstractitensornetwork.jl`):
-  ```julia
-  factorize(tn::AbstractITensorNetwork, edge::AbstractEdge; kws...)
-  factorize(tn::AbstractITensorNetwork, edge::Pair; kws...)
-  ```
-
-* QR-factorize across `edge`, placing the orthogonal factor on `src(edge)` and the remainder
-  on `dst(edge)` (`abstractitensornetwork.jl`):
-  ```julia
-  qr(tn::AbstractITensorNetwork, edge::AbstractEdge; kws...)
-  ```
-
-* SVD-factorize across `edge` (`abstractitensornetwork.jl`):
-  ```julia
-  svd(tn::AbstractITensorNetwork, edge::AbstractEdge; kws...)
+* QR-factorize across `edge`: keep the orthogonal factor on `src(edge)` and absorb
+  the residual into `tn[dst(edge)]`. Vertex set unchanged; the bond between
+  `src(edge)` and `dst(edge)` is re-gauged (`abstractitensornetwork.jl`):
+  ```julia
+  qr(tn::AbstractITensorNetwork, edge::AbstractEdge)
+  qr(tn::AbstractITensorNetwork, edge::Pair)
+  qr!(tn::AbstractITensorNetwork, edge::AbstractEdge)
+  qr!(tn::AbstractITensorNetwork, edge::Pair)
   ```
 
 * Truncate the bond on `edge` via SVD; forwards `cutoff`, `maxdim`, `mindim` kwargs (`abstractitensornetwork.jl`):
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index a0ac6814..9c1e73d5 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -5,10 +5,10 @@ using DataGraphs: DataGraphs, edge_data, get_vertex_data, is_vertex_assigned,
 using Dictionaries: Dictionary
 using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edges, edgetype,
     has_edge, ne, neighbors, rem_edge!, src, vertices
-using ITensors: ITensors, @Algorithm_str, ITensor, addtags, commoninds, commontags,
-    contract, dag, inds, noprime, onehot, prime, replaceprime, replacetags, setprime,
-    settags, sim, swaptags, tags
-using LinearAlgebra: LinearAlgebra, factorize, qr
+using ITensors: ITensors, ITensor, addtags, commoninds, commontags, contract, dag, inds,
+    noprime, onehot, prime, replaceprime, replacetags, setprime, settags, sim, swaptags,
+    tags
+using LinearAlgebra: LinearAlgebra, factorize, qr, qr!
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
@@ -111,14 +111,12 @@ function DataGraphs.get_vertex_data(is::AbstractITensorNetwork, v)
 end
 
 function DataGraphs.set_vertex_data!(tn::AbstractITensorNetwork, value, v)
-    # v = to_vertex(tn, index...)
     @preserve_graph tn[v] = value
     fix_edges!(tn, v)
     return tn
 end
 
 function DataGraphs.set_vertices_data!(tn::AbstractITensorNetwork, values, vertices)
-    # v = to_vertex(tn, index...)
     for v in vertices
         @preserve_graph tn[v] = values[v]
     end
@@ -128,6 +126,13 @@ function DataGraphs.set_vertices_data!(tn::AbstractITensorNetwork, values, verti
     return tn
 end
 
+# Reconcile the graph edges incident to `v` with the shared Indices of
+# `tn[v]`, after a write that may have changed which neighbors share an
+# Index with `v`. Removes any incident edges that no longer correspond
+# to a shared Index, and adds edges to any vertex now sharing one. The
+# long-term invariant is graph-edge ↔ shared-Index one-to-one; until
+# the rest of the codebase maintains it explicitly on every write, this
+# auto-reconciliation runs after `set_vertex_data!` / `set_vertices_data!`.
 function fix_edges!(tn::AbstractITensorNetwork, v)
     for edge in incident_edges(tn, v)
         rem_edge!(tn, edge)
@@ -485,138 +490,36 @@ function ITensors.tags(tn::AbstractITensorNetwork, edge)
     return commontags(is)
 end
 
-function LinearAlgebra.svd(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-    return svd(tn, edgetype(tn)(edge))
-end
-
-function LinearAlgebra.svd(
-        tn::AbstractITensorNetwork,
-        edge::AbstractEdge;
-        U_vertex = src(edge),
-        S_vertex = (edge, "S"),
-        V_vertex = (edge, "V"),
-        u_tags = tags(tn, edge),
-        v_tags = tags(tn, edge),
-        kwargs...
-    )
-    tn = copy(tn)
+# QR-factor `tn[src(edge)]` along the inds it doesn't share with
+# `tn[dst(edge)]`, keep the orthogonal factor `Q` on `src(edge)`, and
+# absorb the residual `R` into `tn[dst(edge)]`. Mutates `tn` in place;
+# the graph is unchanged.
+function LinearAlgebra.qr!(tn::AbstractITensorNetwork, edge::AbstractEdge)
     left_inds = setdiff(inds(tn[src(edge)]), inds(tn[dst(edge)]))
-    U, S, V =
-        svd(tn[src(edge)], left_inds; lefttags = u_tags, righttags = v_tags, kwargs...)
-
-    rem_vertex!(tn, src(edge))
-    add_vertex!(tn, U_vertex)
-    tn[U_vertex] = U
-
-    add_vertex!(tn, S_vertex)
-    tn[S_vertex] = S
-
-    add_vertex!(tn, V_vertex)
-    tn[V_vertex] = V
-
-    return tn
-end
-
-function LinearAlgebra.qr(
-        tn::AbstractITensorNetwork,
-        edge::AbstractEdge;
-        Q_vertex = src(edge),
-        R_vertex = (edge, "R"),
-        tags = tags(tn, edge),
-        kwargs...
+    Q, R = factorize(
+        tn[src(edge)], left_inds; tags = edge_tag(edge), ortho = "left"
     )
-    tn = copy(tn)
-    left_inds = setdiff(inds(tn[src(edge)]), inds(tn[dst(edge)]))
-    Q, R = factorize(tn[src(edge)], left_inds; tags, kwargs...)
-
-    rem_vertex!(tn, src(edge))
-    add_vertex!(tn, Q_vertex)
-    tn[Q_vertex] = Q
-
-    add_vertex!(tn, R_vertex)
-    tn[R_vertex] = R
-
+    @preserve_graph tn[src(edge)] = Q
+    @preserve_graph tn[dst(edge)] = tn[dst(edge)] * R
     return tn
 end
 
-function LinearAlgebra.factorize(
-        tn::AbstractITensorNetwork,
-        edge::AbstractEdge;
-        X_vertex = src(edge),
-        Y_vertex = ("Y", edge),
-        tags = tags(tn, edge),
-        kwargs...
-    )
-    # Promote vertex type
-    V = promote_type(vertextype(tn), typeof(X_vertex), typeof(Y_vertex))
-
-    # TODO: Check `ITensorNetwork{V}`, shouldn't need a copy here.
-    tn = ITensorNetwork{V}(copy(tn))
-
-    neighbors_X = setdiff(neighbors(tn, src(edge)), [dst(edge)])
-    left_inds = setdiff(inds(tn[src(edge)]), inds(tn[dst(edge)]))
-    X, Y = factorize(tn[src(edge)], left_inds; tags, kwargs...)
-
-    rem_vertex!(tn, src(edge))
-    add_vertex!(tn, X_vertex)
-    add_vertex!(tn, Y_vertex)
-
-    add_edge!(data_graph(tn), X_vertex => Y_vertex)
-    for nX in neighbors_X
-        add_edge!(data_graph(tn), X_vertex => nX)
-    end
-    add_edge!(data_graph(tn), Y_vertex => dst(edge))
-    @preserve_graph tn[X_vertex] = X
-    @preserve_graph tn[Y_vertex] = Y
-    return tn
+function LinearAlgebra.qr!(tn::AbstractITensorNetwork, edge::Pair)
+    return qr!(tn, edgetype(tn)(edge))
 end
 
-function LinearAlgebra.factorize(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-    return factorize(tn, edgetype(tn)(edge); kwargs...)
-end
+LinearAlgebra.qr(tn::AbstractITensorNetwork, edge::AbstractEdge) = qr!(copy(tn), edge)
+LinearAlgebra.qr(tn::AbstractITensorNetwork, edge::Pair) = qr!(copy(tn), edge)
 
-# For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function gauge_edge(
-        alg::Algorithm"orthogonalize", tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...
-    )
-    # tn = factorize(tn, edge; kwargs...)
-    # # TODO: Implement as `only(common_neighbors(tn, src(edge), dst(edge)))`
-    # new_vertex = only(neighbors(tn, src(edge)) ∩ neighbors(tn, dst(edge)))
-    # return contract(tn, new_vertex => dst(edge))
-    !has_edge(tn, edge) && throw(ArgumentError("Edge not in graph."))
-    tn = copy(tn)
-    left_inds = setdiff(inds(tn[src(edge)]), inds(tn[dst(edge)]))
-    ltags = tags(tn, edge)
-    X, Y = factorize(tn[src(edge)], left_inds; tags = ltags, ortho = "left", kwargs...)
-    @preserve_graph tn[src(edge)] = X
-    @preserve_graph tn[dst(edge)] = tn[dst(edge)] * Y
-    return tn
-end
-
-# For ambiguity error; TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
-function gauge_walk(
-        alg::Algorithm, tn::AbstractITensorNetwork, edges::Vector{<:AbstractEdge}; kwargs...
-    )
+function gauge_walk(tn::AbstractITensorNetwork, edges)
     tn = copy(tn)
     for edge in edges
-        tn = gauge_edge(alg, tn, edge; kwargs...)
+        qr!(tn, edge)
     end
     return tn
 end
 
-function gauge_walk(alg::Algorithm, tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-    return gauge_edge(alg::Algorithm, tn, edgetype(tn)(edge); kwargs...)
-end
-
-function gauge_walk(
-        alg::Algorithm, tn::AbstractITensorNetwork, edges::Vector{<:Pair}; kwargs...
-    )
-    return gauge_walk(alg, tn, edgetype(tn).(edges); kwargs...)
-end
-
-function tree_gauge(alg::Algorithm, ψ::AbstractITensorNetwork, region)
-    return tree_gauge(alg, ψ, [region])
-end
+tree_gauge(ψ::AbstractITensorNetwork, region) = tree_gauge(ψ, [region])
 
 #Get the path that moves the gauge from a to b on a tree
 #TODO: Move to NamedGraphs
@@ -630,31 +533,18 @@ end
 
 # Gauge a ITensorNetwork from cur_region towards new_region, treating
 # the network as a tree spanned by a spanning tree.
-function tree_gauge(
-        alg::Algorithm,
-        ψ::AbstractITensorNetwork,
-        cur_region::Vector,
-        new_region::Vector;
-        kwargs...
-    )
+function tree_gauge(ψ::AbstractITensorNetwork, cur_region::Vector, new_region::Vector)
     es = edge_sequence_between_regions(ψ, cur_region, new_region)
-    ψ = gauge_walk(alg, ψ, es; kwargs...)
-    return ψ
+    return gauge_walk(ψ, es)
 end
 
 # Gauge a ITensorNetwork towards a region, treating
 # the network as a tree spanned by a spanning tree.
-function tree_gauge(alg::Algorithm, ψ::AbstractITensorNetwork, region::Vector)
-    return tree_gauge(alg, ψ, collect(vertices(ψ)), region)
-end
-
-function tree_orthogonalize(ψ::AbstractITensorNetwork, cur_region, new_region; kwargs...)
-    return tree_gauge(Algorithm("orthogonalize"), ψ, cur_region, new_region; kwargs...)
+function tree_gauge(ψ::AbstractITensorNetwork, region::Vector)
+    return tree_gauge(ψ, collect(vertices(ψ)), region)
 end
 
-function tree_orthogonalize(ψ::AbstractITensorNetwork, region; kwargs...)
-    return tree_gauge(Algorithm("orthogonalize"), ψ, region; kwargs...)
-end
+tree_orthogonalize(ψ::AbstractITensorNetwork, args...) = tree_gauge(ψ, args...)
 
 # TODO: decide whether to use graph mutating methods when resulting graph is unchanged?
 function _truncate_edge(tn::AbstractITensorNetwork, edge::AbstractEdge; kwargs...)
@@ -792,20 +682,6 @@ function linkdims(tn::AbstractITensorNetwork{V}) where {V}
     return ld
 end
 
-# Factor `tn[src(edge)]` via `f` along the inds it doesn't share with
-# `tn[dst(edge)]`, leaving the new factor on `src(edge)` and contracting
-# the residual into `tn[dst(edge)]`. With `f = qr` and no shared inds yet,
-# this threads a fresh link Index between the two endpoints. With shared
-# inds present, it regauges across the existing link.
-function factorize_edge!(f, tn::AbstractITensorNetwork, edge)
-    src_t, dst_t = tn[src(edge)], tn[dst(edge)]
-    left_inds = setdiff(inds(src_t), inds(dst_t))
-    x, y = f(src_t, left_inds)
-    @preserve_graph tn[src(edge)] = x
-    @preserve_graph tn[dst(edge)] = y * dst_t
-    return tn
-end
-
 # Ensure `edge` is present in both the underlying graph and as a shared
 # Index between the endpoint tensors. The long-term invariant is that the
 # two are one-to-one; today they can drift apart, so this fills in
@@ -815,7 +691,7 @@ end
 function Graphs.add_edge!(tn::AbstractITensorNetwork, edge)
     added = !has_edge(tn, edge)
     added && add_edge!(data_graph(tn), edge)
-    isempty(linkinds(tn, edge)) && factorize_edge!(qr, tn, edge)
+    isempty(linkinds(tn, edge)) && qr!(tn, edge)
     return added
 end
 
@@ -854,10 +730,10 @@ function add(tn1::AbstractITensorNetwork, tn2::AbstractITensorNetwork)
     tn1 = copy(tn1)
     tn2 = copy(tn2)
     for e in edges(tn1)
-        length(linkinds(tn1, e)) > 1 && factorize_edge!(qr, tn1, e)
+        length(linkinds(tn1, e)) > 1 && qr!(tn1, e)
     end
     for e in edges(tn2)
-        length(linkinds(tn2, e)) > 1 && factorize_edge!(qr, tn2, e)
+        length(linkinds(tn2, e)) > 1 && qr!(tn2, e)
     end
 
     edges_tn1, edges_tn2 = edges(tn1), edges(tn2)
diff --git a/src/treetensornetworks/abstracttreetensornetwork.jl b/src/treetensornetworks/abstracttreetensornetwork.jl
index 85ebfa18..cf355663 100644
--- a/src/treetensornetworks/abstracttreetensornetwork.jl
+++ b/src/treetensornetworks/abstracttreetensornetwork.jl
@@ -32,17 +32,15 @@ function set_ortho_region(tn::AbstractTTN, new_region)
     return error("Not implemented")
 end
 
-function gauge(alg::Algorithm, ttn::AbstractTTN, region::Vector; kwargs...)
-    ttn = tree_gauge(alg, ttn, collect(ortho_region(ttn)), region; kwargs...)
+function gauge(ttn::AbstractTTN, region::Vector)
+    ttn = tree_gauge(ttn, collect(ortho_region(ttn)), region)
     return set_ortho_region(ttn, region)
 end
 
-function gauge(alg::Algorithm, ttn::AbstractTTN, region; kwargs...)
-    return gauge(alg, ttn, [region]; kwargs...)
-end
+gauge(ttn::AbstractTTN, region) = gauge(ttn, [region])
 
 """
-    orthogonalize(ttn::AbstractTreeTensorNetwork, region; kwargs...) -> TreeTensorNetwork
+    orthogonalize(ttn::AbstractTreeTensorNetwork, region) -> TreeTensorNetwork
 
 Bring `ttn` into an orthogonal gauge with orthogonality center at `region`.
 `region` may be a single vertex or a vector of vertices.
@@ -53,13 +51,9 @@ right-orthogonal with respect to that path.
 
 See also: [`ortho_region`](@ref), [`truncate`](@ref).
 """
-function orthogonalize(ttn::AbstractTTN, region; kwargs...)
-    return gauge(Algorithm("orthogonalize"), ttn, region; kwargs...)
-end
+orthogonalize(ttn::AbstractTTN, region) = gauge(ttn, region)
 
-function tree_orthogonalize(ttn::AbstractTTN, args...; kwargs...)
-    return orthogonalize(ttn, args...; kwargs...)
-end
+tree_orthogonalize(ttn::AbstractTTN, args...) = orthogonalize(ttn, args...)
 
 #
 # Truncation
diff --git a/test/test_itensornetwork.jl b/test/test_itensornetwork.jl
index 05c6dda1..95336570 100644
--- a/test/test_itensornetwork.jl
+++ b/test/test_itensornetwork.jl
@@ -8,7 +8,6 @@ using ITensorNetworks: ITensorNetworks, ITensorNetwork, IndsNetwork, contraction
 using ITensors.NDTensors: NDTensors, dim
 using ITensors: ITensors, ITensor, Index, Op, commonind, commoninds, contract, dag, hasinds,
     inds, inner, itensor, onehot, order, prime, random_itensor, scalartype, sim
-using LinearAlgebra: factorize
 using NDTensors: NDTensors, dim
 using NamedGraphs.GraphsExtensions: disjoint_union, incident_edges
 using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid
@@ -144,12 +143,6 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64})
         rng = StableRNG(1234)
         tn = random_tensornetwork(rng, named_grid(4); link_space = 2)
         Z = norm_sqr(tn)
-        tn_ortho = factorize(tn, 4 => 3)
-        # TODO: Error here in arranging the edges. Arrange by hash?
-        Z̃ = norm_sqr(tn_ortho)
-        @test nv(tn_ortho) == 5
-        @test nv(tn) == 4
-        @test Z ≈ Z̃
         tn_ortho = tree_orthogonalize(tn, [3, 4])
         Z̃ = norm_sqr(tn_ortho)
         @test nv(tn_ortho) == 4

From 47634dfc175dc3e3fde07608d96638c8078c7ea6 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 16:03:22 -0400
Subject: [PATCH 09/11] Delete fix_edges!
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`set_vertex_data!` was auto-reconciling graph edges to match the shared
indices on a written tensor. Inline edge maintenance at the two callsites
that relied on it instead: rewrite the tensors-only `ITensorNetwork{V}`
constructor with explicit O(n²) edge inference after the writes, and
update `_add_edge!` in `test/utils.jl` to use `@preserve_graph` plus an
explicit `add_edge!` on the underlying graph.
---
 src/abstractitensornetwork.jl | 30 ++----------------------------
 src/itensornetwork.jl         | 14 +++++++++++---
 test/utils.jl                 | 13 ++++++++-----
 3 files changed, 21 insertions(+), 36 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 9c1e73d5..2539316d 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -4,7 +4,7 @@ using DataGraphs: DataGraphs, edge_data, get_vertex_data, is_vertex_assigned,
     set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data
 using Dictionaries: Dictionary
 using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edges, edgetype,
-    has_edge, ne, neighbors, rem_edge!, src, vertices
+    has_edge, ne, neighbors, src, vertices
 using ITensors: ITensors, ITensor, addtags, commoninds, commontags, contract, dag, inds,
     noprime, onehot, prime, replaceprime, replacetags, setprime, settags, sim, swaptags,
     tags
@@ -12,7 +12,7 @@ using LinearAlgebra: LinearAlgebra, factorize, qr, qr!
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
-    add_edges, directed_graph, incident_edges, rename_vertices, vertextype, ⊔
+    add_edges, directed_graph, rename_vertices, vertextype, ⊔
 using NamedGraphs: NamedGraphs, NamedGraph, Vertices, not_implemented, steiner_tree
 using SplitApplyCombine: flatten
 
@@ -112,7 +112,6 @@ end
 
 function DataGraphs.set_vertex_data!(tn::AbstractITensorNetwork, value, v)
     @preserve_graph tn[v] = value
-    fix_edges!(tn, v)
     return tn
 end
 
@@ -120,31 +119,6 @@ function DataGraphs.set_vertices_data!(tn::AbstractITensorNetwork, values, verti
     for v in vertices
         @preserve_graph tn[v] = values[v]
     end
-    for v in vertices
-        fix_edges!(tn, v)
-    end
-    return tn
-end
-
-# Reconcile the graph edges incident to `v` with the shared Indices of
-# `tn[v]`, after a write that may have changed which neighbors share an
-# Index with `v`. Removes any incident edges that no longer correspond
-# to a shared Index, and adds edges to any vertex now sharing one. The
-# long-term invariant is graph-edge ↔ shared-Index one-to-one; until
-# the rest of the codebase maintains it explicitly on every write, this
-# auto-reconciliation runs after `set_vertex_data!` / `set_vertices_data!`.
-function fix_edges!(tn::AbstractITensorNetwork, v)
-    for edge in incident_edges(tn, v)
-        rem_edge!(tn, edge)
-    end
-    for vertex in vertices(tn)
-        if v ≠ vertex
-            edge = v => vertex
-            if !isempty(linkinds(tn, edge))
-                add_edge!(data_graph(tn), edge)
-            end
-        end
-    end
     return tn
 end
 
diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl
index e0f6ce38..20c10b9f 100644
--- a/src/itensornetwork.jl
+++ b/src/itensornetwork.jl
@@ -79,8 +79,9 @@ end
 # Construction from collections of ITensors
 #
 
-# Tensors only: derive graph from `keys(tensors)`, then run edge inference.
-# Without the reverse index map, edge inference is O(n²).
+# Tensors only: derive the vertex list from `keys(tensors)`, write the
+# tensors at each vertex, then infer edges from shared Indices in an
+# O(n²) sweep. Without a reverse index map, that's the only available cost.
 function ITensorNetwork{V}(tensors::_ITensorCollection) where {V}
     # Build the vertex list with element type `V` so that an empty `tensors`
     # input doesn't get the graph's vertex type inferred to whatever
@@ -96,7 +97,14 @@ function ITensorNetwork{V}(tensors::_ITensorCollection) where {V}
     default = Dict{V, ITensor}(v => ITensor() for v in vertices(g))
     tn = ITensorNetwork(default, g)
     for v in vertices(g)
-        tn[v] = tensors[v]
+        @preserve_graph tn[v] = tensors[v]
+    end
+    vs = collect(vertices(tn))
+    for i in eachindex(vs), j in (i + 1):lastindex(vs)
+        v1, v2 = vs[i], vs[j]
+        if !isempty(commoninds(tn[v1], tn[v2]))
+            add_edge!(data_graph(tn), v1 => v2)
+        end
     end
     return tn
 end
diff --git a/test/utils.jl b/test/utils.jl
index ffcec865..b86432d9 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -4,9 +4,9 @@
 # inside its gensym module.
 
 using DataGraphs: underlying_graph, vertex_data
-using Graphs: AbstractGraph, dst, edges, src, vertices
+using Graphs: AbstractGraph, add_edge!, dst, edges, src, vertices
 using ITensorNetworks.ITensorsExtensions: trivial_space
-using ITensorNetworks: ITensorNetwork, IndsNetwork
+using ITensorNetworks: @preserve_graph, ITensorNetwork, IndsNetwork, data_graph
 using ITensors.NDTensors: dim
 using ITensors: ITensors, ITensor, Index, dag, itensor, onehot
 using NamedGraphs.GraphsExtensions: incident_edges
@@ -78,12 +78,15 @@ end
 # uses QR, which would push site-state QN flux into the link and leave BP's
 # `default_message` (single-index `delta`) with no compatible block. `onehot`
 # defaults to `Float64`, so we typecast it to `elt` to keep `productstate(elt,
-# ...)` element-type-preserving.
+# ...)` element-type-preserving. Uses `@preserve_graph` plus an explicit
+# `add_edge!` on the underlying graph so the graph-edge ↔ shared-Index
+# invariant is maintained without relying on auto-reconciliation.
 function _add_edge!(elt::Type, tn, edge)
     iₑ = Index(trivial_space(tn), "Link")
     X = ITensors.convert_eltype(elt, onehot(iₑ => 1))
-    tn[src(edge)] *= X
-    tn[dst(edge)] *= dag(X)
+    @preserve_graph tn[src(edge)] = tn[src(edge)] * X
+    @preserve_graph tn[dst(edge)] = tn[dst(edge)] * dag(X)
+    add_edge!(data_graph(tn), edge)
     return tn
 end
 

From 305602940d34ce4968c7fcd8241c09489abd672d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 16:17:34 -0400
Subject: [PATCH 10/11] Duck-type ITensorNetwork ctors, simplify qr!, fix gauge
 docs
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- Drop the `_ITensorCollection` Union from the `ITensorNetwork` constructor
  signatures and the in-file comments referencing the long-removed
  function-callback ctor's dispatch loop (those constructors went away in
  #356). The `Vector{ITensor}` / `Dict{V, ITensor}` / `Dictionary{V, ITensor}`
  shapes still flow through the same body via duck-typed `keys` /
  `keytype` / `tensors[v]`.
- Use `LinearAlgebra.qr` directly inside `qr!` (instead of
  `factorize(...; ortho = "left")`), and reorder the residual contraction
  to `R * tn[dst(edge)]` so the expression reads left-to-right as
  "absorb the residual into the destination tensor". Drop the now-unused
  `factorize` import.
- Refresh the stale `gauge` entry in `docs/src/experimental_methods.md` —
  the `alg::Algorithm` positional and `kwargs...` were dropped from
  `gauge` earlier in this PR.
---
 docs/src/experimental_methods.md |  8 ++++----
 src/abstractitensornetwork.jl    |  8 +++-----
 src/itensornetwork.jl            | 26 ++++----------------------
 3 files changed, 11 insertions(+), 31 deletions(-)

diff --git a/docs/src/experimental_methods.md b/docs/src/experimental_methods.md
index de602acc..cd8a5abd 100644
--- a/docs/src/experimental_methods.md
+++ b/docs/src/experimental_methods.md
@@ -53,10 +53,10 @@ Methods which still need to be discussed, modified, or deprecated.
   tree-traversal QRs; `orthogonalize` is the user-facing wrapper, with `tree_orthogonalize`
   as an alias (`treetensornetworks/abstracttreetensornetwork.jl`):
   ```julia
-  gauge(alg::Algorithm, ttn::AbstractTTN, region::Vector; kwargs...)
-  gauge(alg::Algorithm, ttn::AbstractTTN, region; kwargs...)
-  orthogonalize(ttn::AbstractTTN, region; kwargs...)
-  tree_orthogonalize(ttn::AbstractTTN, args...; kwargs...)
+  gauge(ttn::AbstractTTN, region::Vector)
+  gauge(ttn::AbstractTTN, region)
+  orthogonalize(ttn::AbstractTTN, region)
+  tree_orthogonalize(ttn::AbstractTTN, args...)
   ```
 
 * Sweep-based truncation. The whole-TTN form orthogonalizes towards `src(e)` before
diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 2539316d..706ab8a2 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -8,7 +8,7 @@ using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edge
 using ITensors: ITensors, ITensor, addtags, commoninds, commontags, contract, dag, inds,
     noprime, onehot, prime, replaceprime, replacetags, setprime, settags, sim, swaptags,
     tags
-using LinearAlgebra: LinearAlgebra, factorize, qr, qr!
+using LinearAlgebra: LinearAlgebra, qr, qr!
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
@@ -470,11 +470,9 @@ end
 # the graph is unchanged.
 function LinearAlgebra.qr!(tn::AbstractITensorNetwork, edge::AbstractEdge)
     left_inds = setdiff(inds(tn[src(edge)]), inds(tn[dst(edge)]))
-    Q, R = factorize(
-        tn[src(edge)], left_inds; tags = edge_tag(edge), ortho = "left"
-    )
+    Q, R = qr(tn[src(edge)], left_inds; tags = edge_tag(edge))
     @preserve_graph tn[src(edge)] = Q
-    @preserve_graph tn[dst(edge)] = tn[dst(edge)] * R
+    @preserve_graph tn[dst(edge)] = R * tn[dst(edge)]
     return tn
 end
 
diff --git a/src/itensornetwork.jl b/src/itensornetwork.jl
index 20c10b9f..defa489c 100644
--- a/src/itensornetwork.jl
+++ b/src/itensornetwork.jl
@@ -40,21 +40,13 @@ julia> tn = ITensorNetwork([ITensor(i, j), ITensor(j, k)]);
 
 See also: `IndsNetwork`, [`TreeTensorNetwork`](@ref ITensorNetworks.TreeTensorNetwork).
 """
-const _ITensorCollection = Union{
-    AbstractVector{<:ITensor},
-    AbstractDict{<:Any, <:ITensor},
-    AbstractDictionary{<:Any, <:ITensor},
-}
-
 struct ITensorNetwork{V} <: AbstractITensorNetwork{V}
     data_graph::DataGraph{V, ITensor, ITensor, NamedGraph{V}, NamedEdge{V}}
 
     # Sole inner ctor: place `tensors` at the vertices of `graph`. No checks —
     # `tensors` must be indexable at every vertex, the graph's edges are
     # taken at face value.
-    function ITensorNetwork{V}(
-            tensors::_ITensorCollection, graph::NamedGraph
-        ) where {V}
+    function ITensorNetwork{V}(tensors, graph::NamedGraph) where {V}
         g = NamedGraph{V}(graph)
         dg = DataGraph(g; vertex_data_type = ITensor, edge_data_type = ITensor)
         for v in vertices(g)
@@ -82,18 +74,11 @@ end
 # Tensors only: derive the vertex list from `keys(tensors)`, write the
 # tensors at each vertex, then infer edges from shared Indices in an
 # O(n²) sweep. Without a reverse index map, that's the only available cost.
-function ITensorNetwork{V}(tensors::_ITensorCollection) where {V}
+function ITensorNetwork{V}(tensors) where {V}
     # Build the vertex list with element type `V` so that an empty `tensors`
     # input doesn't get the graph's vertex type inferred to whatever
     # `keys(tensors)` happens to give (e.g. `Int` for an empty `Vector{ITensor}`).
     g = NamedGraph(V[v for v in keys(tensors)])
-    # Annotate the default Dict's element type explicitly: when the
-    # comprehension is empty (e.g. constructing an empty `ITensorNetwork{V}`
-    # via the bootstrap call inside the IndsNetwork function-callback ctor)
-    # type inference can give `Dict{Any, Any}` for `V = Any`, which then
-    # falls outside `_ITensorCollection` and routes the next dispatch
-    # through the wrong path (function-callback Group C) — infinite
-    # recursion. Locking the value type to `ITensor` avoids that.
     default = Dict{V, ITensor}(v => ITensor() for v in vertices(g))
     tn = ITensorNetwork(default, g)
     for v in vertices(g)
@@ -110,10 +95,10 @@ function ITensorNetwork{V}(tensors::_ITensorCollection) where {V}
 end
 
 # Non-parametric delegates: extract `V` via `keytype` / `vertextype`.
-function ITensorNetwork(tensors::_ITensorCollection)
+function ITensorNetwork(tensors)
     return ITensorNetwork{keytype(tensors)}(tensors)
 end
-function ITensorNetwork(tensors::_ITensorCollection, graph::NamedGraph)
+function ITensorNetwork(tensors, graph::NamedGraph)
     return ITensorNetwork{vertextype(graph)}(tensors, graph)
 end
 
@@ -123,9 +108,6 @@ end
 
 function ITensorNetwork{V}(tn::ITensorNetwork) where {V}
     g = NamedGraph{V}(underlying_graph(tn))
-    # Type-annotated so empty-`tn` + `V = Any` doesn't drop us out of
-    # `_ITensorCollection` and into the wrong dispatch (see the parametric
-    # `(tensors)` form above for the same gotcha).
     tensors = Dict{V, ITensor}(v => tn[v] for v in vertices(tn))
     return ITensorNetwork(tensors, g)
 end

From a9205b45a7c79b2b41af6bf5d9e97919b1197869 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 11 May 2026 16:29:30 -0400
Subject: [PATCH 11/11] Restore fix_edges! as the set_*_data!
 auto-reconciliation hook
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Plain `tn[v] = ...` writes go through `set_vertex_data!`, which needs
to keep the graph edges incident to `v` in sync with the shared indices
of the new tensor. Restore `fix_edges!` for that purpose, keyed off the
non-`@preserve_graph` path; the two callsites that now maintain the
invariant explicitly (the tensors-only `ITensorNetwork{V}` constructor
and `_add_edge!` in `test/utils.jl`) bypass `set_vertex_data!` via
`@preserve_graph` and don't pay the cost. Plan is to drop this helper
once the invariant becomes structural — e.g., once an `ITensorNetwork`
reverse Index map makes auto-reconciliation either unnecessary or a
constant-time lookup.
---
 src/abstractitensornetwork.jl | 33 +++++++++++++++++++++++++++++++--
 1 file changed, 31 insertions(+), 2 deletions(-)

diff --git a/src/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index 706ab8a2..3b51701c 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -4,7 +4,7 @@ using DataGraphs: DataGraphs, edge_data, get_vertex_data, is_vertex_assigned,
     set_vertex_data!, underlying_graph, underlying_graph_type, vertex_data
 using Dictionaries: Dictionary
 using Graphs: Graphs, Graph, add_edge!, add_vertex!, bfs_tree, center, dst, edges, edgetype,
-    has_edge, ne, neighbors, src, vertices
+    has_edge, ne, neighbors, rem_edge!, src, vertices
 using ITensors: ITensors, ITensor, addtags, commoninds, commontags, contract, dag, inds,
     noprime, onehot, prime, replaceprime, replacetags, setprime, settags, sim, swaptags,
     tags
@@ -12,7 +12,7 @@ using LinearAlgebra: LinearAlgebra, qr, qr!
 using MacroTools: @capture
 using NDTensors: NDTensors, Algorithm, dim, scalartype
 using NamedGraphs.GraphsExtensions:
-    add_edges, directed_graph, rename_vertices, vertextype, ⊔
+    add_edges, directed_graph, incident_edges, rename_vertices, vertextype, ⊔
 using NamedGraphs: NamedGraphs, NamedGraph, Vertices, not_implemented, steiner_tree
 using SplitApplyCombine: flatten
 
@@ -112,6 +112,7 @@ end
 
 function DataGraphs.set_vertex_data!(tn::AbstractITensorNetwork, value, v)
     @preserve_graph tn[v] = value
+    fix_edges!(tn, v)
     return tn
 end
 
@@ -119,6 +120,34 @@ function DataGraphs.set_vertices_data!(tn::AbstractITensorNetwork, values, verti
     for v in vertices
         @preserve_graph tn[v] = values[v]
     end
+    for v in vertices
+        fix_edges!(tn, v)
+    end
+    return tn
+end
+
+# Reconcile the graph edges incident to `v` with the shared Indices of
+# `tn[v]`, after a non-`@preserve_graph` write that may have changed
+# which neighbors share an Index with `v`. Removes any incident edges
+# that no longer correspond to a shared Index, and adds edges to any
+# vertex now sharing one. The long-term invariant is graph-edge ↔
+# shared-Index one-to-one; until that's structural (e.g. via a reverse
+# Index map on the type), this auto-reconciliation runs after
+# `set_vertex_data!` / `set_vertices_data!` so plain `tn[v] = ...`
+# writes don't desync the two. Callers that already maintain the
+# invariant explicitly should bypass this with `@preserve_graph`.
+function fix_edges!(tn::AbstractITensorNetwork, v)
+    for edge in incident_edges(tn, v)
+        rem_edge!(tn, edge)
+    end
+    for vertex in vertices(tn)
+        if v ≠ vertex
+            edge = v => vertex
+            if !isempty(linkinds(tn, edge))
+                add_edge!(data_graph(tn), edge)
+            end
+        end
+    end
     return tn
 end