diff --git a/docs/src/experimental_methods.md b/docs/src/experimental_methods.md
index 0b7ca2ad..cd8a5abd 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
@@ -62,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
@@ -224,12 +215,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/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/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
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/abstractitensornetwork.jl b/src/abstractitensornetwork.jl
index ed98ed45..3b51701c 100644
--- a/src/abstractitensornetwork.jl
+++ b/src/abstractitensornetwork.jl
@@ -5,14 +5,14 @@ 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 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, qr, 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
 
@@ -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,16 @@ 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 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)
@@ -485,138 +493,34 @@ 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
-
+    Q, R = qr(tn[src(edge)], left_inds; tags = edge_tag(edge))
+    @preserve_graph tn[src(edge)] = Q
+    @preserve_graph tn[dst(edge)] = R * tn[dst(edge)]
     return tn
 end
 
-function LinearAlgebra.qr(
-        tn::AbstractITensorNetwork,
-        edge::AbstractEdge;
-        Q_vertex = src(edge),
-        R_vertex = (edge, "R"),
-        tags = tags(tn, edge),
-        kwargs...
-    )
-    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
-
-    return tn
+function LinearAlgebra.qr!(tn::AbstractITensorNetwork, edge::Pair)
+    return qr!(tn, edgetype(tn)(edge))
 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
-end
-
-function LinearAlgebra.factorize(tn::AbstractITensorNetwork, edge::Pair; kwargs...)
-    return factorize(tn, edgetype(tn)(edge); kwargs...)
-end
-
-# 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
+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_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 +534,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)
+function tree_gauge(ψ::AbstractITensorNetwork, region::Vector)
+    return tree_gauge(ψ, collect(vertices(ψ)), region)
 end
 
-function tree_orthogonalize(ψ::AbstractITensorNetwork, cur_region, new_region; kwargs...)
-    return tree_gauge(Algorithm("orthogonalize"), ψ, cur_region, new_region; kwargs...)
-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...)
@@ -709,38 +600,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
@@ -824,40 +683,17 @@ 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,
-# 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
+# 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)) && qr!(tn, edge)
+    return added
 end
 
 # TODO: What to output? Could be an `IndsNetwork`. Or maybe
@@ -874,10 +710,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
 
@@ -894,15 +726,23 @@ 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 && qr!(tn1, e)
+    end
+    for e in edges(tn2)
+        length(linkinds(tn2, e)) > 1 && qr!(tn2, e)
+    end
 
     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)
+        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/itensornetwork.jl b/src/itensornetwork.jl
index e0f6ce38..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)
@@ -79,33 +71,34 @@ 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²).
-function ITensorNetwork{V}(tensors::_ITensorCollection) where {V}
+# 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) 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)
-        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
 
 # 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
 
@@ -115,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
diff --git a/src/tebd.jl b/src/tebd.jl
index 5f4891f1..a9de4d00 100644
--- a/src/tebd.jl
+++ b/src/tebd.jl
@@ -15,11 +15,11 @@ function tebd(
     s = siteinds(ψ)
     u⃗ = Vector{ITensor}(𝒰, s)
     nsteps = Int(β ÷ Δβ)
+    ψ = add_edges(ψ, edges(ψ))
     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(ψ)
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/src/treetensornetworks/treetensornetwork.jl b/src/treetensornetworks/treetensornetwork.jl
index 7dc7cb25..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,46 +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)
-    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)]))
-        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)
-        a = a_r
-    end
-    tn[ortho_center] = a
-    ttn_a = TreeTensorNetwork(tn)
-    return orthogonalize(ttn_a, ortho_center)
-end
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
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)
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