From 7f9ee757ecf3ec8ed53b73a2a8c9eef7a5ce82bb Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Sun, 5 Jul 2026 15:47:53 -0400
Subject: [PATCH 01/14] Use graded-capable matrix functions for gauging,
 entanglement, and the compat eigen

Simple update gauges environments with `sqrth_invsqrth_safe` and `renyi_entropy` uses `sqrth_safe`, replacing the eigendecomposition helpers vendored in `utils.jl`. The compat `eigen` projects onto the Hermitian part with `MatrixAlgebraKit.project_hermitian` at the matricized level instead of a named-level add. The compat factorizations bipartition indices by name so graded bond duality does not drop legs. `Rxxyy` and `Rxxyyzz` are built as whole matrices from raising and lowering operators so they project onto symmetric backends, and the vendored gate set gains the `SWAP` matrix. Removes the singular value sign fix after gate application (MatrixAlgebraKit singular values are nonnegative, so it was a no-op). Fermionic sign handling is left for future work. Dense `expect` on networks with pruned trivial tensors is fixed by pruning at the network level. Adds `TensorAlgebra` as a direct dependency.

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
---
 Project.toml                                  |   8 +-
 src/Apply/apply_gates.jl                      |  12 +-
 src/Apply/gate_definitions.jl                 |  16 ++-
 src/Apply/simple_update.jl                    |   7 +-
 .../ITensorsITensorBaseCompat.jl              |   2 +-
 src/ITensorsITensorBaseCompat/itensors.jl     | 114 ++++++++++--------
 src/ITensorsITensorBaseCompat/ops.jl          |  74 ++++++++----
 src/MessagePassing/beliefpropagationcache.jl  |   5 +-
 src/MessagePassing/loopcorrection.jl          |   5 +-
 src/TensorNetworks/abstracttensornetwork.jl   |  26 +---
 src/TensorNetworks/tensornetworkstate.jl      |  32 ++++-
 src/contraction_sequences.jl                  |  14 ++-
 src/entanglement.jl                           |   9 +-
 src/imports.jl                                |   6 +-
 src/utils.jl                                  |  26 +---
 test/test_contraction_sequences.jl            |  10 +-
 16 files changed, 209 insertions(+), 157 deletions(-)

diff --git a/Project.toml b/Project.toml
index 672c0ad..1aa3886 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,6 +21,7 @@ SimpleGraphConverter = "205b04f2-f585-4877-a239-566270b3f673"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 
@@ -28,7 +29,7 @@ TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138"
 Adapt = "4.3.0"
 Dictionaries = "0.4"
 Graphs = "1.8.0"
-ITensorBase = "0.10.2"
+ITensorBase = "0.10.9"
 KrylovKit = "0.10.2"
 LinearAlgebra = "1.11.0"
 MatrixAlgebraKit = "0.6.8"
@@ -40,6 +41,7 @@ SimpleGraphConverter = "0.1.0"
 SplitApplyCombine = "1.2.3"
 Statistics = "1.11.1"
 StatsBase = "0.34.4"
+TensorAlgebra = "0.17"
 TensorOperations = "5.2"
 TypeParameterAccessors = "0.3.10, 0.4"
 julia = "1.10"
@@ -50,3 +52,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "Random"]
+
+[sources]
+ITensorBase = {path = "/Users/mfishman/.julia/dev/ITensorBase/.worktrees/mf/tensoralgebra-0.17"}
+TensorAlgebra = {path = "/Users/mfishman/.julia/dev/TensorAlgebra/.worktrees/mf/project-trailing-axes"}
diff --git a/src/Apply/apply_gates.jl b/src/Apply/apply_gates.jl
index 2225e71..0f0d1f8 100644
--- a/src/Apply/apply_gates.jl
+++ b/src/Apply/apply_gates.jl
@@ -111,11 +111,13 @@ function apply_gate!(
     if length(v⃗) == 2
         v1, v2 = v⃗
         e = NamedEdge(v1 => v2)
-        ind2 = commonind(s_values, first(updated_tensors))
-        δuv = dag(copy(s_values))
-        δuv = replaceind(δuv, ind2, prime(ind2))
-        map_diag!(sign, δuv, δuv)
-        s_values = denseblocks(s_values) * denseblocks(δuv)
+        u = commonind(s_values, first(updated_tensors))
+        v = noncommonind(s_values, first(updated_tensors))
+        # The new messages are the singular values over the bond-and-prime pair.
+        # MatrixAlgebraKit singular values are nonnegative, so the legacy sign fix
+        # (`s * sign(s)` via `map_diag!`) was a no-op and is dropped; fermionic sign
+        # handling for this message construction is future work.
+        s_values = replaceind(s_values, v, prime(u))
         setmessage!(ψ_bpc, e, dag(s_values))
         setmessage!(ψ_bpc, reverse(e), s_values)
     end
diff --git a/src/Apply/gate_definitions.jl b/src/Apply/gate_definitions.jl
index 7762de4..c348fcd 100644
--- a/src/Apply/gate_definitions.jl
+++ b/src/Apply/gate_definitions.jl
@@ -263,8 +263,13 @@ ITensors.op(o::OpName"xx_plus_yy", ::SiteType"Qubit"; θ::Number, β::Number) =
 
 Gate for rotation by XXYY at a given angle.
 """
-function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2", s1::Index, s2::Index; θ = 1)
-    h = 0.5 * (op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2))
+function ITensors.op(::OpName"Rxxyy", ::SiteType"S=1/2"; θ = 1)
+    # Built as one two-site matrix in the manifestly charge-conserving σ± form,
+    # ½(XX + YY) = σ⁺σ⁻ + σ⁻σ⁺, rather than from single-site `op("X", s)` factors:
+    # the gate conserves U(1) charge, but a standalone `X` does not, so the factored
+    # construction has no symmetric representation even though the sum does.
+    σp, σm = [0.0 1.0; 0.0 0.0], [0.0 0.0; 1.0 0.0]
+    h = kron(σp, σm) + kron(σm, σp)
     return exp(-im * θ * h)
 end
 ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) =
@@ -275,8 +280,11 @@ ITensors.op(o::OpName"Rxxyy", ::SiteType"Qubit"; θ::Number) =
 
 Gate for rotation by XXYYZZ at a given angle.
 """
-function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2", s1::Index, s2::Index; θ = 1)
-    h = 0.5 * (op("X", s1) * op("X", s2) + op("Y", s1) * op("Y", s2) + op("Z", s1) * op("Z", s2))
+function ITensors.op(::OpName"Rxxyyzz", ::SiteType"S=1/2"; θ = 1)
+    # One two-site matrix in the σ± form, for the same reason as `Rxxyy` above:
+    # ½(XX + YY + ZZ) = σ⁺σ⁻ + σ⁻σ⁺ + ½ ZZ.
+    σp, σm, σz = [0.0 1.0; 0.0 0.0], [0.0 0.0; 1.0 0.0], [1.0 0.0; 0.0 -1.0]
+    h = kron(σp, σm) + kron(σm, σp) + 0.5 * kron(σz, σz)
     return exp(-im * θ * h)
 end
 ITensors.op(o::OpName"Rxxyyzz", ::SiteType"Qubit"; θ::Number) =
diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl
index 40b7644..d7aa292 100644
--- a/src/Apply/simple_update.jl
+++ b/src/Apply/simple_update.jl
@@ -35,8 +35,11 @@ function simple_update(
         envs_v2 = filter(env -> hascommoninds(env, ψ⃗[2]), envs)
         @assert all(ndims(env) == 2 for env in vcat(envs_v1, envs_v2))
 
-        sqrt_inv_sqrt_envs_v1 = pseudo_sqrt_inv_sqrt.(envs_v1; cutoff = sqrt_cutoff)
-        sqrt_inv_sqrt_envs_v2 = pseudo_sqrt_inv_sqrt.(envs_v2; cutoff = sqrt_cutoff)
+        sqrt_invsqrt = env -> sqrth_invsqrth_safe(
+            env, (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0
+        )
+        sqrt_inv_sqrt_envs_v1 = map(sqrt_invsqrt, envs_v1)
+        sqrt_inv_sqrt_envs_v2 = map(sqrt_invsqrt, envs_v2)
         sqrt_envs_v1, inv_sqrt_envs_v1 = first.(sqrt_inv_sqrt_envs_v1), last.(sqrt_inv_sqrt_envs_v1)
         sqrt_envs_v2, inv_sqrt_envs_v2 = first.(sqrt_inv_sqrt_envs_v2), last.(sqrt_inv_sqrt_envs_v2)
 
diff --git a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
index da3b3b5..b32f293 100644
--- a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
+++ b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
@@ -30,7 +30,7 @@ export
     # Index operations
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
     # ITensor construction
-    itensor, random_itensor, scalar, delta, onehot, combiner, combinedind,
+    itensor, random_itensor, scalar, delta, onehot,
     # Factorizations
     qr, svd, eigen, factorize, factorize_svd,
     # Diagonal manipulation
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index d9200b0..4540681 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -3,9 +3,11 @@
 #
 # Strategy (see ITensorDevelopmentPlans api_migration_map.md):
 #   - Names below are thin wrappers over ITensorBase / TensorAlgebra / MatrixAlgebraKit.
-#   - `combiner`, the factorization return shapes, `map_diag`, the operator/SiteType
-#     system, and the boundary-MPS (ITensorMPS) paths are NOT wrapped here; they need
-#     callsite translation or upstream stack work and are tracked separately.
+#   - The factorization return shapes, `map_diag`, the operator/SiteType system, and
+#     the boundary-MPS (ITensorMPS) paths are NOT wrapped here; they need callsite
+#     translation or upstream stack work and are tracked separately. The legacy
+#     `combiner` is retired rather than wrapped: call sites use the next-gen fusion
+#     primitives (`matricize` below, the fused identity `Base.one`) directly.
 #
 # ITensorBase keeps most of this API internal (unexported), so we reach for the
 # qualified names and re-publish the legacy spellings into the TNQS namespace.
@@ -15,6 +17,7 @@ using Adapt: Adapt
 using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange, name,
     nameddims, plev, tags, unnamed
 using LinearAlgebra: LinearAlgebra
+using TensorAlgebra.MatrixAlgebra: MatrixAlgebra
 
 # Legacy `inds(t; plev, tags)` accepted index-filtering keywords; the next-gen
 # `ITensorBase.inds` takes none. Own a compat `inds` that forwards to `ITensorBase.inds`
@@ -300,10 +303,11 @@ function _hermitian_eigh(
     )
     cod, dom = _astuple(codomain), _astuple(domain)
     # `MAK.eigh_full` rejects a matrix that is hermitian only up to numerical noise, but
-    # the caller asserts `ishermitian = true`. Project onto the hermitian part first
-    # (the conjugate transpose swaps codomain ↔ domain), matching legacy `eigen`'s
-    # treat-as-hermitian behavior.
-    m = (m + replaceinds(conj(m), (cod..., dom...), (dom..., cod...))) / 2
+    # the caller asserts `ishermitian = true`. Project onto the hermitian part first,
+    # matching legacy `eigen`'s treat-as-hermitian behavior. `MAK.project_hermitian`
+    # works at the matricized level, so it stays generic over graded and `TensorMap`
+    # backends (a named-level `m + swap(conj(m))` is not).
+    m = MAK.project_hermitian(m, cod, dom)
     D, U = MAK.eigh_full(m, cod, dom)
     u = only(commoninds(D, U))         # eigenvalue index shared with U
     t = only(uniqueinds(D, U))         # D's other (independent) index
@@ -312,8 +316,8 @@ function _hermitian_eigh(
 end
 
 # Partitioned form `eigen(m, Linds, Rinds)` reproduces legacy ITensors' reconstruction
-# `m = Vt * D * dag(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`), which
-# is what `eigendecomp`/`pseudo_sqrt_inv_sqrt` rely on — no conjugation of `U`.
+# `m = Vt * D * dag(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`) —
+# no conjugation of `U`.
 function eigen(m::AbstractITensor, codomain, domain; kwargs...)
     return _hermitian_eigh(m, codomain, domain; kwargs...)
 end
@@ -349,6 +353,19 @@ function _trunc_spec(cutoff, maxdim)
     return reduce(&, specs)
 end
 
+# Split `a`'s indices into (left, right) by the requested `linds`, matching by name:
+# the caller's `Index` objects may carry the other endpoint's (dual) space on a graded
+# backend, so `Index`-equality filtering silently drops them (see the name-matching
+# note on the index-set algebra above). Both groups are `a`'s own indices, so the
+# spaces handed to the factorization are `a`'s.
+function _bipartition_inds(a::AbstractITensor, linds)
+    lnames = name.(cat_inds(linds...))
+    allinds = collect(inds(a))
+    left = filter(i -> name(i) ∈ lnames, allinds)
+    right = filter(i -> name(i) ∉ lnames, allinds)
+    return left, right
+end
+
 # Legacy `factorize(a, linds...; ortho, cutoff, maxdim, tags)` splits `a` into `L * R`
 # with the new bond between them. `linds` selects `L`'s indices (besides the bond) and
 # may be passed splatted, as a single index, or as one collection. `ortho = "left"`
@@ -363,9 +380,7 @@ function factorize(
         maxdim = nothing,
         tags = nothing
     )
-    allinds = collect(inds(a))
-    left = filter(∈(allinds), cat_inds(linds...))
-    right = setdiff(allinds, left)
+    left, right = _bipartition_inds(a, linds)
     trunc = _trunc_spec(cutoff, maxdim)
     if isnothing(trunc)
         if ortho == "left"
@@ -412,8 +427,10 @@ function _absorb_svd(U, S, Vt, ortho)
         u = only(commoninds(U, S))
         v = only(commoninds(S, Vt))
         up = prime(u)
-        sqrtσ = sqrt.(abs.(diagview(unnamed(S))))
-        return U * diagonaltensor(sqrtσ, (u, up)), Vt * diagonaltensor(sqrtσ, (v, up))
+        # `S` is diagonal, so this hits the diagonal fast path (no eigendecomposition)
+        # on every backend. No clamping: zero singular values stay zero under `^(1//2)`.
+        sqrtS = MatrixAlgebra.sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0)
+        return U * replaceind(sqrtS, v, up), Vt * replaceind(sqrtS, u, up)
     end
     return error(
         "compat `factorize_svd` supports ortho = \"left\" / \"right\" / \"none\" (got $(repr(ortho)))"
@@ -433,9 +450,7 @@ function factorize_svd(
         maxdim = nothing,
         tags = nothing
     )
-    allinds = collect(inds(a))
-    left = filter(∈(allinds), cat_inds(linds...))
-    right = setdiff(allinds, left)
+    left, right = _bipartition_inds(a, linds)
     trunc = _trunc_spec(cutoff, maxdim)
     U, S, Vt = if isnothing(trunc)
         MAK.svd_compact(a, Tuple(left), Tuple(right))
@@ -443,7 +458,7 @@ function factorize_svd(
         MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc)
     end
     total = abs2(LinearAlgebra.norm(a))
-    kept = sum(abs2, diagview(unnamed(S)))
+    kept = abs2(LinearAlgebra.norm(S))
     truncerr = if iszero(total)
         zero(real(scalartype(a)))
     else
@@ -460,53 +475,46 @@ function factorize_svd(
 end
 
 #
-# Combiner / index fusion. Legacy `combiner(is...)` returns a tensor that, when
-# contracted, fuses `is` into one index; `dag(C)` splits it back, and
-# `combinedind(C)` recovers that fused index. Reproduced densely as a reshaped
-# identity (`C[is..., c] = δ(c, flatten(is))`), which the next-gen fusion
-# primitives back via ordinary contraction. The fused index carries a sentinel tag
-# so `combinedind` can find it. Dense only; a graded combiner is a stack gap.
-const _COMBINER_TAG = "combiner"
-function combiner(is::Tuple; eltype::Type = Float64)
-    d = prod(length, is)
-    c = ITensorBase.settag(Index(d), _COMBINER_TAG, "combined")
-    return reshape(Matrix{eltype}(LinearAlgebra.I, d, d), (length.(is)..., d))[is..., c]
-end
-combiner(is::Index...; kwargs...) = combiner(is; kwargs...)
-combiner(is::AbstractVector{<:Index}; kwargs...) = combiner(Tuple(is); kwargs...)
-function combinedind(c::AbstractITensor)
-    return only(
-        filter(i -> get(ITensorBase.tags(i), _COMBINER_TAG, "") == "combined", inds(c))
-    )
-end
+# Index fusion. The legacy `combiner` is retired (no compat shim); call sites fuse
+# index groups with the next-gen `matricize(t, row_inds => row_name, col_inds =>
+# col_name)` (minting each fused name via `name(uniquename(i))`) or build a fused
+# identity with `Base.one`, both of which are graded-capable. `matricize` is the
+# `TensorAlgebra` generic, extended by ITensorBase for named tensors.
+using TensorAlgebra: matricize
 
 #
 # Diagonal manipulation. Legacy `map_diag(f, T)` applies `f` to the diagonal of a
 # (diagonal-like) tensor; used on factorization spectra (singular values /
-# eigenvalues). Reproduced via the same diagonal machinery as `delta`.
-_diagcartesian(arr, k) = CartesianIndex(ntuple(Returns(k), ndims(arr)))
+# eigenvalues). Reproduced via `MAK.diagview`, which aliases the diagonal storage on
+# every backend (a plain view for dense, blockwise for graded, per-sector for a
+# `TensorMap`), so no scalar indexing is needed.
+function _map_diagview!(f, dv)
+    copyto!(dv, map(f, dv))
+    return dv
+end
+# A backend's `diagview` may be a dict of per-block diagonal views (e.g. a plain
+# `TensorMap`, keyed by sector) rather than a single vector view.
+function _map_diagview!(f, dv::AbstractDict)
+    foreach(v -> _map_diagview!(f, v), values(dv))
+    return dv
+end
 function map_diag(f, T::AbstractITensor)
     arr = copy(unnamed(T))
-    for k in 1:minimum(size(arr))
-        idx = _diagcartesian(arr, k)
-        arr[idx] = f(arr[idx])
-    end
-    return arr[inds(T)...]
+    _map_diagview!(f, MAK.diagview(arr))
+    return nameddims(arr, ITensorBase.dimnames(T))
 end
 function map_diag!(f, T::AbstractITensor)
-    arr = unnamed(T)
-    for k in 1:minimum(size(arr))
-        idx = _diagcartesian(arr, k)
-        arr[idx] = f(arr[idx])
-    end
+    _map_diagview!(f, MAK.diagview(unnamed(T)))
     return T
 end
 # Out-of-place-into-`dest` form `map_diag!(f, dest, src)`: write `f` of `src`'s diagonal
 # onto `dest`'s diagonal (TNQS calls it with `dest === src` for an in-place diagonal map).
 function map_diag!(f, dest::AbstractITensor, src::AbstractITensor)
     d, s = unnamed(dest), unnamed(src)
-    for k in 1:minimum(size(s))
-        d[_diagcartesian(d, k)] = f(s[_diagcartesian(s, k)])
+    if d === s
+        map_diag!(f, dest)
+    else
+        copyto!(MAK.diagview(d), map(f, MAK.diagview(s)))
     end
     return dest
 end
@@ -532,6 +540,10 @@ data(T::AbstractITensor) = unnamed(T)
 # alone, so reproduce the eltype conversion for a `Number` target (TNQS uses
 # `adapt(eltype)(state(...))` to build typed product states).
 function Adapt.adapt_structure(::Type{elt}, T::AbstractITensor) where {elt <: Number}
+    # Short-circuit the no-op case: a non-`AbstractArray` backend (e.g. a `TensorMap`)
+    # has no `convert(AbstractArray{elt}, ...)` method, but needs none when the element
+    # type already matches.
+    eltype(T) === elt && return T
     return nameddims(convert(AbstractArray{elt}, unnamed(T)), ITensorBase.dimnames(T))
 end
 
diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/ITensorsITensorBaseCompat/ops.jl
index 079564c..c5f35a0 100644
--- a/src/ITensorsITensorBaseCompat/ops.jl
+++ b/src/ITensorsITensorBaseCompat/ops.jl
@@ -10,6 +10,8 @@
 # only ever calls `op(name, sites...; kwargs...)` and `state(name, site)`, so the vendor
 # is a plain name-keyed lookup rather than a reimplementation of that dispatch machinery.
 
+using TensorAlgebra: TensorAlgebra, project, tryproject
+
 #
 # Named single-site states. `state(name, i)` returns the named state vector as an
 # `ITensor` over `i`.
@@ -20,15 +22,28 @@ const _STATE_VECTORS = Dict{String, Vector{ComplexF64}}(
     "X+" => [1, 1] / sqrt(2), "+" => [1, 1] / sqrt(2),
     "X-" => [1, -1] / sqrt(2), "-" => [1, -1] / sqrt(2),
     "Y+" => [1, im] / sqrt(2),
-    "Y-" => [1, -im] / sqrt(2),
+    "Y-" => [1, -im] / sqrt(2)
 )
 function state(name::AbstractString, i::Index)
     v = get(_STATE_VECTORS, name) do
-        error("unknown single-site state \"$name\" (vendored states: $(sort(collect(keys(_STATE_VECTORS)))))")
+        return error(
+            "unknown single-site state \"$name\" (vendored states: $(sort(collect(keys(_STATE_VECTORS)))))"
+        )
     end
     length(v) == length(i) ||
-        error("state \"$name\" has dimension $(length(v)) but the site index has dimension $(length(i))")
-    return collect(v)[i]
+        error(
+        "state \"$name\" has dimension $(length(v)) but the site index has dimension $(length(i))"
+    )
+    ψ = tryproject(v, (i,))
+    isnothing(ψ) || return ψ
+    # The state carries a charge under the site index's grading (e.g. "Dn" on a
+    # U(1) site), so it can't live in the flux-zero space over `i` alone. Carry the
+    # charge on an explicit length-1 auxiliary leg: project with a trailing axis so
+    # the backend derives the leg's sector from the state vector, then wrap the
+    # derived axis in a freshly named `Index`.
+    raw = project(reshape(v, (length(v), 1)), (unnamed(i),), ())
+    aux = Index(TensorAlgebra.axes(raw, 2))
+    return nameddims(raw, (ITensorBase.name(i), ITensorBase.name(aux)))
 end
 
 #
@@ -53,13 +68,15 @@ macro SiteType_str(name)
 end
 
 # Embed a `d^n × d^n` operator matrix (computational basis, first site most
-# significant) into an `ITensor` over `(prime.(sites)..., sites...)`.
+# significant) into an `ITensor` with codomain `prime.(sites)` (outputs) and
+# domain `sites` (inputs). `project` is checked, so on graded sites an operator
+# that is not symmetric under the site index's grading throws an `InexactError`.
 function _op_matrix_to_itensor(M::AbstractMatrix, sites::Tuple)
     ds = length.(sites)
     n = length(sites)
     A = reshape(Matrix{ComplexF64}(M), (reverse(ds)..., reverse(ds)...))
     A = permutedims(A, (reverse(1:n)..., reverse((n + 1):(2n))...))
-    return A[ITensorBase.prime.(sites)..., sites...]
+    return project(A, ITensorBase.prime.(sites), sites)
 end
 
 # Top-level `op(name, sites...; kwargs...)`. The identity is dimension-general (used
@@ -68,7 +85,10 @@ end
 function op(name::AbstractString, sites::Index...; kwargs...)
     if name in ("I", "Id") && length(sites) == 1
         s = only(sites)
-        return Matrix{ComplexF64}(LinearAlgebra.I, length(s), length(s))[ITensorBase.prime(s), s]
+        return project(
+            Matrix{ComplexF64}(LinearAlgebra.I, length(s), length(s)),
+            (ITensorBase.prime(s),), (s,)
+        )
     end
     return op(OpName(name), SiteType("S=1/2"), sites...; kwargs...)
 end
@@ -87,15 +107,17 @@ const _GATE_MATRICES = Dict{String, Matrix{ComplexF64}}(
     "Z" => [1 0; 0 -1],
     "H" => [1 1; 1 -1] / sqrt(2),
     "S" => [1 0; 0 im],
-    "T" => [1 0; 0 cis(π / 4)],
+    "T" => [1 0; 0 cis(π / 4)]
 )
 function _gate_matrix(name::AbstractString; kwargs...)
     name in ("Rx", "Ry", "Rz", "P") || return get(_GATE_MATRICES, name) do
-        error("unknown single-site operator \"$name\" (vendored operators: $(sort(collect(keys(_GATE_MATRICES)))) plus Rx/Ry/Rz/P)")
+        return error(
+            "unknown single-site operator \"$name\" (vendored operators: $(sort(collect(keys(_GATE_MATRICES)))) plus Rx/Ry/Rz/P)"
+        )
     end
     if name == "Rx"
         θ = kwargs[:θ]
-        return ComplexF64[cos(θ / 2) (-im*sin(θ / 2)); (-im*sin(θ / 2)) cos(θ / 2)]
+        return ComplexF64[cos(θ / 2) (-im * sin(θ / 2)); (-im * sin(θ / 2)) cos(θ / 2)]
     elseif name == "Ry"
         θ = kwargs[:θ]
         return ComplexF64[cos(θ / 2) (-sin(θ / 2)); sin(θ / 2) cos(θ / 2)]
@@ -108,7 +130,9 @@ function _gate_matrix(name::AbstractString; kwargs...)
     end
 end
 # Matrix method for any single-qubit name handled by `_gate_matrix`.
-op(::OpName{Name}, ::SiteType"S=1/2"; kwargs...) where {Name} = _gate_matrix(String(Name); kwargs...)
+function op(::OpName{Name}, ::SiteType"S=1/2"; kwargs...) where {Name}
+    return _gate_matrix(String(Name); kwargs...)
+end
 
 # Two-qubit gate matrices (computational basis, first site most significant).
 const _σx = ComplexF64[0 1; 1 0]
@@ -117,25 +141,23 @@ const _σz = ComplexF64[1 0; 0 -1]
 op(::OpName"Rxx", ::SiteType"S=1/2"; ϕ) = exp(-im * ϕ * kron(_σx, _σx))
 op(::OpName"Ryy", ::SiteType"S=1/2"; ϕ) = exp(-im * ϕ * kron(_σy, _σy))
 op(::OpName"Rzz", ::SiteType"S=1/2"; ϕ) = exp(-im * ϕ * kron(_σz, _σz))
-op(::OpName"CPHASE", ::SiteType"S=1/2"; ϕ) = ComplexF64[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 cis(ϕ)]
+function op(::OpName"CPHASE", ::SiteType"S=1/2"; ϕ)
+    return ComplexF64[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 cis(ϕ)]
+end
+function op(::OpName"SWAP", ::SiteType"S=1/2")
+    return Float64[1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
+end
 
-# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor` (matricize over the
-# index/prime pairs, exponentiate, rebuild). Legacy ITensors provided `exp(::ITensor)`; gates
-# defined as `exp` of a Hamiltonian operator rely on it (e.g. a user
+# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor`, inferring the
+# prime-pair codomain/domain and forwarding to ITensorBase's matricization
+# `exp(a, dimnames_codomain, dimnames_domain)` (graded-capable). Legacy ITensors provided
+# `exp(::ITensor)`; gates defined as `exp` of a Hamiltonian operator rely on it (e.g. a user
 # `op(::OpName"MyZRot", ...) = exp(-im θ/2 * op("Z", s))`). To de-pirate: make this compat-owned
-# (an `exp` in this module, not a `Base.exp` method), inferring the prime-pair codomain/domain and
-# forwarding to ITensorBase's matricization `exp(a, dimnames_codomain, dimnames_domain)` — which
-# already exists and is graded-capable, so this dense combiner-based version goes away. Not an
-# upstream candidate (the upstream matricization `exp` is the target, not a `Base.exp(::ITensor)`).
+# (an `exp` in this module, not a `Base.exp` method). Not an upstream candidate (the upstream
+# matricization `exp` is the target, not a `Base.exp(::ITensor)`).
 function Base.exp(t::AbstractITensor)
     p0 = filter(i -> ITensorBase.plev(i) == 0, collect(inds(t)))
     isempty(p0) && error("exp(::ITensor) expects indices paired as (i, prime(i))")
     p1 = map(ITensorBase.prime, p0)
-    cr, cc = combiner(Tuple(p1)), combiner(Tuple(p0))
-    tc = (t * cr) * cc
-    rci, cci = combinedind(cr), combinedind(cc)
-    d = length(rci)
-    M = ComplexF64[tc[rci => a, cci => b] for a in 1:d, b in 1:d]
-    E = exp(M)[rci, cci]
-    return (E * dag(cr)) * dag(cc)
+    return exp(t, Tuple(p1), Tuple(p0))
 end
diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl
index f864504..93fbc70 100644
--- a/src/MessagePassing/beliefpropagationcache.jl
+++ b/src/MessagePassing/beliefpropagationcache.jl
@@ -170,8 +170,9 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge}
     seq = contraction_sequence(tensors; alg = "omeinsum", optimizer = GreedyMethod())
     t = contract(tensors; sequence = seq)
 
-    row_combiner, col_combiner = ITensors.combiner(e_virtualinds), ITensors.combiner(e_virtualinds_sim)
-    t = t * row_combiner * col_combiner
+    row_name = name(ITensorBase.uniquename(first(e_virtualinds)))
+    col_name = name(ITensorBase.uniquename(first(e_virtualinds_sim)))
+    t = ITensors.matricize(t, Tuple(e_virtualinds) => row_name, Tuple(e_virtualinds_sim) => col_name)
     t = adapt(Vector{ComplexF64})(t)
     t = ITensors.array(t)
     λs = reverse(sort(LinearAlgebra.eigvals(t); by = abs))
diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl
index cc5298c..a37a59c 100644
--- a/src/MessagePassing/loopcorrection.jl
+++ b/src/MessagePassing/loopcorrection.jl
@@ -50,10 +50,7 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg)
                     row_inds = vcat(row_inds, dag.(prime.(row_inds)))
                     col_inds = vcat(col_inds, dag.(prime.(col_inds)))
                 end
-                row_combiner, col_combiner = combiner(row_inds), combiner(col_inds)
-                ap =
-                    adapt_like(message(bpc, e), denseblocks(delta(combinedind(col_combiner), dag(combinedind(row_combiner)))))
-                ap = ap * row_combiner * dag(col_combiner)
+                ap = adapt_like(message(bpc, e), identity_tensor(row_inds, col_inds))
                 ap = ap - message(bpc, e) * mer
                 push!(antiprojectors, ap)
             end
diff --git a/src/TensorNetworks/abstracttensornetwork.jl b/src/TensorNetworks/abstracttensornetwork.jl
index 0203abb..444f66f 100644
--- a/src/TensorNetworks/abstracttensornetwork.jl
+++ b/src/TensorNetworks/abstracttensornetwork.jl
@@ -28,13 +28,15 @@ function maxvirtualdim(tn::AbstractTensorNetwork)
     return maximum(maximum.([dim.(virtualinds(tn, e)) for e in edges(tn)]))
 end
 
+# Compare by name, not by `Index` equality: a shared graded link is stored nondual
+# on one endpoint and dual (conjugated) on the other, so the two `Index` objects
+# differ even though they name the same bond.
 function ITensors.uniqueinds(tn::AbstractTensorNetwork, v)
     tv_inds = Index[i for i in inds(tn[v])]
     vns = neighbors(tn, v)
     isempty(vns) && return tv_inds
-    neighbor_inds = reduce(vcat, [Index[i for i in inds(tn[vn])] for vn in vns])
-    is = setdiff(tv_inds, neighbor_inds)
-    return is
+    neighbor_names = reduce(vcat, [[name(i) for i in inds(tn[vn])] for vn in vns])
+    return filter(i -> name(i) ∉ neighbor_names, tv_inds)
 end
 
 function setindex_preserve!(tn::AbstractTensorNetwork, value::ITensor, vertex)
@@ -105,24 +107,6 @@ function map_virtualinds(f::Function, tn::AbstractTensorNetwork)
     return map_virtualinds!(f, tn)
 end
 
-function combine_virtualinds!(tn::AbstractTensorNetwork)
-    dtype = datatype(tn)
-    for e in edges(tn)
-        vinds = ITensors.commoninds(tn[src(e)], tn[dst(e)])
-        if length(vinds) > 1
-            C = adapt(dtype)(ITensors.combiner(vinds))
-            setindex_preserve!(tn, tn[src(e)] * C, src(e))
-            setindex_preserve!(tn, tn[dst(e)] * C, dst(e))
-        end
-    end
-    return tn
-end
-
-function combine_virtualinds(tn::AbstractTensorNetwork)
-    tn = copy(tn)
-    return combine_virtualinds!(tn)
-end
-
 """Add two tensornetworks together. The network structures need to be have the same graph structure"""
 function add(tn1::AbstractTensorNetwork, tn2::AbstractTensorNetwork)
     @assert graph(tn1) == graph(tn2)
diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl
index b4cbc44..fe5d260 100644
--- a/src/TensorNetworks/tensornetworkstate.jl
+++ b/src/TensorNetworks/tensornetworkstate.jl
@@ -22,6 +22,14 @@ tensors(tns::TensorNetworkState) = tensors(tensornetwork(tns))
 Base.copy(tns::TensorNetworkState) = TensorNetworkState(copy(tensornetwork(tns)), copy(siteinds(tns)))
 
 TensorNetworkState(tn::TensorNetwork) = TensorNetworkState(tn, siteinds(tn))
+# Widen a concretely-typed site-index dictionary to the field type.
+function TensorNetworkState(tn::TensorNetwork{V}, sinds::Dictionary) where {V}
+    s = Dictionary{V, Vector{<:Index}}()
+    for v in keys(sinds)
+        set!(s, v, sinds[v])
+    end
+    return TensorNetworkState(tn, s)
+end
 TensorNetworkState(tensors::Dictionary, g::NamedGraph) = TensorNetworkState(TensorNetwork(tensors, g))
 TensorNetworkState(tensors::Union{Dictionary, Vector{<:ITensor}}) = TensorNetworkState(TensorNetwork(tensors))
 
@@ -53,6 +61,12 @@ function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Functi
         sinds = siteinds(tns, v)
         tnv = tns[v]
         tnv_dag = dag(prime(tnv))
+        # Dangling non-physical legs (e.g. a charged state's charge leg) always
+        # pair between ket and bra, so unprime them on the bra unconditionally.
+        auxinds = Index[i for i in setdiff(uniqueinds(tns, v), sinds)]
+        if !isempty(auxinds)
+            tnv_dag = replaceinds(tnv_dag, Index[prime(i) for i in auxinds], auxinds)
+        end
         if op_strings(v) == "ρ" || isempty(sinds)
             append!(factors, ITensor[tnv, tnv_dag])
         elseif op_strings(v) == "I"
@@ -69,9 +83,13 @@ end
 norm_factors(tns::TensorNetworkState, v; kwargs...) = norm_factors(tns, [v]; kwargs...)
 bp_factors(tns::TensorNetworkState, v) = norm_factors(tns, v)
 
+# The flat starting message is the identity between the ket links and their bra
+# copies, built as an identity operator so it follows the links' backend (graded
+# links give a graded message; the legacy `delta` filled a dense diagonal).
 function default_message(tns::TensorNetworkState, edge::AbstractEdge)
     linds = virtualinds(tns, edge)
-    return adapt_like(tns, denseblocks(delta(vcat(linds, prime(dag(linds))))))
+    cod, dom = Tuple(linds), Tuple(prime(dag(linds)))
+    return adapt_like(tns, one(zeros(scalartype(tns), cod..., dom...), cod, dom))
 end
 
 """
@@ -152,12 +170,16 @@ function tensornetworkstate(eltype, f::Function, g::AbstractGraph, siteinds::Dic
         end
     end
 
-    l = Dict(e => Index(1) for e in edges(g))
+    # Trivial links, minted to match the site indices' backend (`trivialrange` of a
+    # graded range is the length-1 trivial-sector range) so graded product states
+    # stay graded. The dst side takes the conjugate so the pair contracts to 1.
     for e in edges(g)
-        tensors[src(e)] *= onehot(eltype, l[e] => 1)
-        tensors[dst(e)] *= onehot(eltype, l[e] => 1)
+        l = Index(trivialrange(unnamed(only(siteinds[src(e)]))))
+        x = fill!(similar(tensors[src(e)], (l,)), true)
+        tensors[src(e)] *= x
+        tensors[dst(e)] *= conj(x)
     end
-    return TensorNetworkState(tensors, g)
+    return TensorNetworkState(TensorNetwork(tensors, g), siteinds)
 end
 
 """
diff --git a/src/contraction_sequences.jl b/src/contraction_sequences.jl
index f3f5732..bcbe3c6 100644
--- a/src/contraction_sequences.jl
+++ b/src/contraction_sequences.jl
@@ -10,9 +10,13 @@ using OMEinsumContractionOrders: OMEinsumContractionOrders, optimize_code, EinCo
 # untouched originals), so no placeholder tensor is needed.
 is_trivial_tensor(t::ITensor) = all(i -> dim(i) == 1, inds(t))
 
+# The sequence optimizers only use indices as opaque labels (plus their dimension), so
+# hand them the index names: a shared leg is stored nondual on one tensor and dual on the
+# other, and comparing whole `Index` objects would see two different labels (or, for a
+# space-backed index, fail to compare at all).
 function contraction_network(tensors::Vector{<:ITensor}; prune_tensors = false)
     return map(tensors) do t
-        is = collect(inds(t))
+        is = collect(name.(inds(t)))
         (prune_tensors && is_trivial_tensor(t)) ? empty(is) : is
     end
 end
@@ -20,7 +24,7 @@ end
 function contraction_sequence(::Algorithm"optimal", tensors::Vector{<:ITensor}; prune_tensors = false)
     network = contraction_network(tensors; prune_tensors)
     #Converting dims to Float64 to minimize overflow issues
-    inds_to_dims = Dict(i => Float64(dim(i)) for i in unique(Iterators.flatten(network)))
+    inds_to_dims = Dict(name(i) => Float64(dim(i)) for t in tensors for i in inds(t))
     seq, _ = optimaltree(network, inds_to_dims)
     seq = typeof(seq) <: Int ? [seq] : seq
     return seq
@@ -38,10 +42,10 @@ end
 
 #OMEinsumContractionOrders helpers
 function to_eincode(tensors::Vector{<:ITensor})
-    ixs = map(t -> collect(inds(t)), tensors)
+    ixs = map(t -> collect(name.(inds(t))), tensors)
     LT = eltype(eltype(ixs))
-    iy = collect(LT, reduce(noncommoninds, tensors))
-    size_dict = Dict{LT, Int}(i => dim(i) for ix in ixs for i in ix)
+    iy = collect(LT, name.(reduce(noncommoninds, tensors)))
+    size_dict = Dict{LT, Int}(name(i) => dim(i) for t in tensors for i in inds(t))
     return EinCode(ixs, iy), size_dict
 end
 
diff --git a/src/entanglement.jl b/src/entanglement.jl
index 33c0f34..c766271 100644
--- a/src/entanglement.jl
+++ b/src/entanglement.jl
@@ -30,7 +30,9 @@ end
 
 function matricize(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a)))
     col_inds = prime.(row_inds)
-    return ITensors.array(a * ITensors.combiner(row_inds) * ITensors.combiner(col_inds))
+    row_name = name(ITensorBase.uniquename(first(row_inds)))
+    col_name = name(ITensorBase.uniquename(first(col_inds)))
+    return ITensors.array(ITensors.matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name))
 end
 
 """
@@ -78,7 +80,10 @@ function renyi_entropy(
     ee = 0
     m1, m2 = message(bp_cache, e), message(bp_cache, reverse(e))
     edge_ind = only(virtualinds(bp_cache, e))
-    root_m2 = first(pseudo_sqrt_inv_sqrt(m2))
+    root_m2 = sqrth_safe(
+        m2, (inds(m2)[1],), (inds(m2)[2],);
+        atol = 10 * eps(real(scalartype(m2))), rtol = 0
+    )
 
     edge_ind_p, edge_ind_pp = prime(edge_ind), prime(prime(edge_ind))
     ρ = (m1 * replaceind(root_m2, edge_ind_p, edge_ind_pp)) * root_m2
diff --git a/src/imports.jl b/src/imports.jl
index a804b80..f140faf 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -56,7 +56,7 @@ import .ITensorsITensorBaseCompat as ITensors
 using .ITensorsITensorBaseCompat:
     inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds,
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
-    itensor, random_itensor, scalar, delta, onehot, combiner, combinedind,
+    itensor, random_itensor, scalar, delta, onehot,
     qr, svd, eigen, factorize, factorize_svd,
     map_diag, map_diag!,
     scalartype, datatype, array, data,
@@ -66,7 +66,9 @@ using .ITensorsITensorBaseCompat:
     Algorithm, @Algorithm_str,
     hastags,
     state, op, OpName, SiteType, @OpName_str, @SiteType_str
-using ITensorBase: ITensorBase, Index, ITensor, name, plev, tags
+using ITensorBase: ITensorBase, Index, ITensor, name, plev, tags, unnamed
+using TensorAlgebra: trivialrange
+using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe
 
 using Adapt: adapt
 
diff --git a/src/utils.jl b/src/utils.jl
index 3e61306..1d39e0e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -15,32 +15,14 @@ function is_ring_graph(g::AbstractGraph)
     return is_line_graph(g_mod)
 end
 
-function pseudo_sqrt_inv_sqrt(M::ITensor; cutoff = 10 * eps(real(scalartype(M))))
-    @assert length(inds(M)) == 2
-    Q, D, Qdag = eigendecomp(M, inds(M)[1], inds(M)[2]; ishermitian = true)
-    D_sqrt = ITensors.map_diag(x -> iszero(x) || abs(x) < cutoff ? 0 : sqrt(x), D)
-    D_inv_sqrt = ITensors.map_diag(x -> iszero(x) || abs(x) < cutoff ? 0 : inv(sqrt(x)), D)
-    M_sqrt = Q * D_sqrt * Qdag
-    M_inv_sqrt = Q * D_inv_sqrt * Qdag
-    return M_sqrt, M_inv_sqrt
-end
-
-#TODO: Make this work for non-hermitian A
-function eigendecomp(A::ITensor, linds, rinds; ishermitian = false, kwargs...)
-    @assert ishermitian
-    D, U = safe_eigen(A, linds, rinds; ishermitian, kwargs...)
-    ul, ur = noncommonind(D, U), commonind(D, U)
-    Ul = replaceinds(U, vcat([rinds], [ur]), vcat([linds], [ul]))
-    return Ul, D, dag(U)
-end
-
 # Adapt `t` to the storage datatype (eltype + device) of `ref`.
 adapt_like(ref, t) = adapt(datatype(ref))(t)
 
+# The fused identity between the two index groups, via `Base.one` (graded-capable);
+# `one` takes a prototype tensor for the names/eltype/spaces, which `zeros` supplies.
 function identity_tensor(eltype, row_inds::Vector{<:Index}, col_inds::Vector{<:Index})
-    c_row, c_col = ITensors.combiner(row_inds),ITensors.combiner(col_inds)
-    t= ITensors.denseblocks(ITensors.delta(eltype, ITensors.combinedind(c_row), ITensors.combinedind(c_col)))
-    return (t * c_row)*c_col
+    row_is, col_is = Tuple(row_inds), Tuple(col_inds)
+    return one(zeros(eltype, row_is, col_is), row_is, col_is)
 end
 
 identity_tensor(row_inds::Vector{<:Index}, col_inds::Vector{<:Index}) = identity_tensor(Float64, row_inds, col_inds)
diff --git a/test/test_contraction_sequences.jl b/test/test_contraction_sequences.jl
index 0ff76f6..da0e44b 100644
--- a/test/test_contraction_sequences.jl
+++ b/test/test_contraction_sequences.jl
@@ -1,5 +1,5 @@
 @eval module $(gensym())
-using ITensorBase: Index
+using ITensorBase: Index, name
 using Random
 using TensorNetworkQuantumSimulator
 const TNQS = TensorNetworkQuantumSimulator
@@ -18,13 +18,15 @@ collect_leaves!(acc, x) = (for y in x; collect_leaves!(acc, y); end; acc)
 
     # --- to_eincode: ITensors -> (EinCode, size_dict). Tests the omeinsum-specific
     #     input conversion directly, so a silent fallback to another backend can't pass it.
+    #     Labels are index names: a shared leg's `Index` differs between its two tensors
+    #     under a graded backend (nondual vs dual), so names are the backend-stable label.
     i, j, k = Index(2), Index(3), Index(4)
     A = ITensors.random_itensor(i, j)
     B = ITensors.random_itensor(j, k)
     code, size_dict = TNQS.to_eincode([A, B])
-    @test Set(Set.(getixsv(code))) == Set([Set([i, j]), Set([j, k])])  # per-tensor index sets
-    @test Set(getiyv(code)) == Set([i, k])                             # open indices (j is contracted)
-    @test size_dict == Dict(i => 2, j => 3, k => 4)
+    @test Set(Set.(getixsv(code))) == Set([Set(name.([i, j])), Set(name.([j, k]))])  # per-tensor index sets
+    @test Set(getiyv(code)) == Set(name.([i, k]))                    # open indices (j is contracted)
+    @test size_dict == Dict(name(i) => 2, name(j) => 3, name(k) => 4)
 
     # --- to_contraction_sequence: NestedEinsum -> nested tensor-position tree. Tests our
     #     converter on hand-built trees with known shapes (deterministic, exact).

From 2b21e2634addcf1379dd540ade50e87e5c474023 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Sun, 5 Jul 2026 17:13:35 -0400
Subject: [PATCH 02/14] Support symmetric backends in inner, rdm, expect, and
 sampling

The order-2 compat `delta` is now built with a checked `project` and carries exactly the indices it is given, matching the legacy convention where callers set the arrow directions. Dangling non-physical legs (for example the charge leg of a charged product state) pair directly between ket and bra, handled by a shared `bra_tensor` helper used by `norm_factors`, the forms, and the samplers. The compat `onehot` carries its charge on a derived auxiliary leg when needed, the compat `array` densifies any backend storage, and the boundary MPS message stitching mints its virtual legs as trivial sectors of the messages' backend.
---
 src/Forms/bilinearform.jl                 |  5 ++-
 src/Forms/quadraticform.jl                |  2 +-
 src/ITensorsITensorBaseCompat/itensors.jl | 42 +++++++++++++++++++----
 src/MessagePassing/boundarympscache.jl    | 11 ++++--
 src/TensorNetworks/tensornetworkstate.jl  | 20 +++++++----
 src/imports.jl                            |  2 +-
 src/rdm.jl                                |  3 +-
 src/sampling.jl                           | 18 ++++++----
 8 files changed, 74 insertions(+), 29 deletions(-)

diff --git a/src/Forms/bilinearform.jl b/src/Forms/bilinearform.jl
index 9dfbdc1..7c397d9 100644
--- a/src/Forms/bilinearform.jl
+++ b/src/Forms/bilinearform.jl
@@ -14,12 +14,11 @@ Base.copy(blf::BilinearForm) = BilinearForm(copy(blf.ket), copy(blf.operator), c
 
 #Constructor, bra is taken to be in the vector space of ket so the dual is taken
 function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState)
-    dtype = datatype(ket)
     @assert graph(ket) == graph(bra)
-    bra = map_tensors(t -> dag(prime(t)), bra)
     sinds = siteinds(ket)
     verts = collect(vertices(ket))
-    operator_tensors = [adapt(dtype)(reduce(*, ITensor[denseblocks(delta(sind, prime(dag(sind)))) for sind in sinds[v]])) for v in verts]
+    bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts]))
+    operator_tensors = [reduce(*, ITensor[delta(scalartype(ket), dag(sind), prime(sind)) for sind in sinds[v]]) for v in verts]
     operator = TensorNetworkState(Dictionary(verts, operator_tensors))
     return BilinearForm(ket, operator, bra)
 end
diff --git a/src/Forms/quadraticform.jl b/src/Forms/quadraticform.jl
index 12f4b0e..5c6ab3a 100644
--- a/src/Forms/quadraticform.jl
+++ b/src/Forms/quadraticform.jl
@@ -6,7 +6,7 @@ end
 ket(qf::QuadraticForm) = qf.ket
 operator(qf::QuadraticForm) = qf.operator
 bra(qf::QuadraticForm) = prime(dag(ket(qf)))
-bra_tensor(qf::QuadraticForm, v) = dag(prime(ket(qf)[v]))
+bra_tensor(qf::QuadraticForm, v) = bra_tensor(ket(qf), v)
 bra_virtualinds(qf::QuadraticForm, edge::NamedEdge) = dag.(prime.(virtualinds(ket(qf), edge)))
 
 Base.copy(qf::QuadraticForm) = QuadraticForm(copy(qf.ket), copy(qf.operator))
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 4540681..933be68 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -18,6 +18,7 @@ using ITensorBase: ITensorBase, AbstractITensor, ITensor, Index, NamedUnitRange,
     nameddims, plev, tags, unnamed
 using LinearAlgebra: LinearAlgebra
 using TensorAlgebra.MatrixAlgebra: MatrixAlgebra
+using TensorAlgebra: TensorAlgebra
 
 # Legacy `inds(t; plev, tags)` accepted index-filtering keywords; the next-gen
 # `ITensorBase.inds` takes none. Own a compat `inds` that forwards to `ITensorBase.inds`
@@ -217,6 +218,20 @@ function diagonaltensor(
 end
 
 delta(eltype::Type, is::Tuple) = diagonaltensor(ones(eltype, minimum(length, is)), is)
+# The order-2 delta is the identity map over the index pair, built via checked
+# `project` so the index axes select the backend (dense, graded, `TensorMap`). The
+# result carries exactly the given indices with their own axes (`project` dualizes
+# domain axes, so `dag(j)` in the domain slot passes `j`'s axis through) — the legacy
+# QN-`delta` convention, where callers arrow the indices themselves. The pair must be
+# arrow-opposite (axis of `j` the dual of axis of `i`, e.g. `(dag(i), prime(i))`); an
+# identity between incompatible gradings throws. Higher-order deltas stay on the dense
+# `diagonaltensor` path above (a super-diagonal generally cannot be embedded while
+# preserving a nontrivial symmetry).
+function delta(eltype::Type, is::Tuple{Index, Index})
+    i, j = is
+    id = Matrix{eltype}(LinearAlgebra.I, length(i), length(j))
+    return TensorAlgebra.project(id, (i,), (dag(j),))
+end
 delta(eltype::Type, is::Index...) = delta(eltype, is)
 delta(eltype::Type, is::AbstractVector{<:Index}) = delta(eltype, Tuple(is))
 delta(is::Tuple) = delta(Float64, is)
@@ -232,16 +247,25 @@ delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is))
 function tr(t::AbstractITensor)
     out = copy(t)
     for i in inds(t; plev = 0)
-        out *= Adapt.adapt(datatype(t), delta(i, prime(i)))
+        # The delta legs carry the duals of `t`'s prime pair (`i` and its primed dual),
+        # so each leg contracts its partner on any backend.
+        out *= delta(scalartype(t), dag(i), prime(i))
     end
     return scalar(out)
 end
 
-# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`).
+# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), built via checked
+# `project` so it follows the index's backend. When the position's sector is charged
+# under a graded `i` the vector can't live over `i` alone, so carry the charge on a
+# derived length-1 auxiliary leg — the same construction as a charged `state` (`ops.jl`).
 function onehot(eltype::Type, (i, p)::Pair{<:Index})
     v = zeros(eltype, length(i))
     v[p] = one(eltype)
-    return v[i]
+    ψ = TensorAlgebra.tryproject(v, (i,))
+    isnothing(ψ) || return ψ
+    raw = TensorAlgebra.project(reshape(v, (length(v), 1)), (unnamed(i),), ())
+    aux = Index(TensorAlgebra.axes(raw, 2))
+    return nameddims(raw, (ITensorBase.name(i), ITensorBase.name(aux)))
 end
 onehot(p::Pair{<:Index}) = onehot(Float64, p)
 
@@ -523,9 +547,12 @@ end
 # Storage / element type accessors.
 #
 # `scalartype` is re-exported above. `datatype` is the underlying storage array type
-# (used by `adapt`); `array` / `data` expose the plain unnamed array.
+# (used by `adapt`); `data` exposes the plain unnamed array. `array` densifies (legacy
+# `array` materialized a dense array from any storage): a no-op on a dense backend,
+# while graded / `TensorMap` storage converts through its canonical flat basis, so
+# positions agree with `onehot` / `project` on the same axes.
 datatype(T::AbstractITensor) = typeof(unnamed(T))
-array(T::AbstractITensor) = unnamed(T)
+array(T::AbstractITensor) = convert(Array, unnamed(T))
 data(T::AbstractITensor) = unnamed(T)
 
 # TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): extends
@@ -648,7 +675,10 @@ end
 # argument). Kept here for now; the call sites get modernized to the `tags` kwarg later, retiring
 # these methods rather than upstreaming them.
 #
-# Legacy positional tagged-index constructor `Index(dim, "tag")`.
+# Legacy positional tagged-index constructor `Index(dim, "tag")`. Only the dimension
+# form is pirated: a range/space first argument collides with ITensorBase's own
+# two-argument `Index` constructors, so tagging an index minted over a backend axis
+# spells the two steps (`settags(Index(r), "tag")`).
 ITensorBase.Index(dim::Integer, tagstr::AbstractString) = settags(Index(dim), tagstr)
 # Build a fresh index carrying a tag dictionary (legacy `Index(dim, tags(i))`, where
 # the next-gen `tags` returns a `Dict{String, String}`).
diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl
index 47e09c8..17ea521 100644
--- a/src/MessagePassing/boundarympscache.jl
+++ b/src/MessagePassing/boundarympscache.jl
@@ -191,11 +191,16 @@ function set_interpartition_messages!(
         end
         for i in 1:(length(es) - 1)
             virt_dim = virtual_index_dimension(bmps_cache, es[i], es[i + 1])
-            ind = Index(virt_dim, "m$(i)$(i + 1)")
             m1, m2 = message(bmps_cache, es[i]), message(bmps_cache, es[i + 1])
-            t = adapt_like(m1, dense(delta(ind)))
+            # The stitching leg is minted trivial (all weight in the charge-0 sector)
+            # so it follows the messages' backend; the all-ones filling matches the
+            # legacy dense `delta(ind)`.
+            ind = ITensors.settags(Index(trivialrange(unnamed(first(inds(m1))), virt_dim)), "m$(i)$(i + 1)")
+            t = fill!(similar(m1, (ind,)), true)
+            # The two copies of the stitching leg contract against each other along the
+            # message MPS, so one side takes the conjugate.
             setmessage!(bmps_cache, es[i], m1 * t)
-            setmessage!(bmps_cache, es[i + 1], m2 * t)
+            setmessage!(bmps_cache, es[i + 1], m2 * dag(t))
         end
     end
     return bmps_cache
diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl
index fe5d260..75cbf83 100644
--- a/src/TensorNetworks/tensornetworkstate.jl
+++ b/src/TensorNetworks/tensornetworkstate.jl
@@ -55,18 +55,24 @@ function Base.setindex!(tns::TensorNetworkState, value::ITensor, v)
     return tns
 end
 
+# Bra copy of a tensor: `dag` and prime all legs except `auxinds` — dangling
+# non-physical legs (e.g. a charged state's charge leg), which always pair directly
+# between ket and bra rather than through an operator or a message.
+function bra_tensor(t::ITensor, auxinds::Vector{<:Index})
+    tdag = dag(prime(t))
+    return isempty(auxinds) ? tdag : replaceinds(tdag, prime.(auxinds), auxinds)
+end
+bra_tensor(tns::TensorNetworkState, v) = bra_tensor(tns[v], auxinds(tns, v))
+
+# The dangling non-physical legs of a vertex tensor: dangling legs that are not site indices.
+auxinds(tns::TensorNetworkState, v) = Index[i for i in setdiff(uniqueinds(tns, v), siteinds(tns, v))]
+
 function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I")
     factors = ITensor[]
     for v in verts
         sinds = siteinds(tns, v)
         tnv = tns[v]
-        tnv_dag = dag(prime(tnv))
-        # Dangling non-physical legs (e.g. a charged state's charge leg) always
-        # pair between ket and bra, so unprime them on the bra unconditionally.
-        auxinds = Index[i for i in setdiff(uniqueinds(tns, v), sinds)]
-        if !isempty(auxinds)
-            tnv_dag = replaceinds(tnv_dag, Index[prime(i) for i in auxinds], auxinds)
-        end
+        tnv_dag = bra_tensor(tns, v)
         if op_strings(v) == "ρ" || isempty(sinds)
             append!(factors, ITensor[tnv, tnv_dag])
         elseif op_strings(v) == "I"
diff --git a/src/imports.jl b/src/imports.jl
index f140faf..5a58c24 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -54,7 +54,7 @@ using TensorOperations
 # referenced `ITensors`, so the per-file imports mirror the original ITensors-based code.
 import .ITensorsITensorBaseCompat as ITensors
 using .ITensorsITensorBaseCompat:
-    inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds,
+    inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, cat_inds,
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
     itensor, random_itensor, scalar, delta, onehot,
     qr, svd, eigen, factorize, factorize_svd,
diff --git a/src/rdm.jl b/src/rdm.jl
index bf942e0..f07c6c4 100644
--- a/src/rdm.jl
+++ b/src/rdm.jl
@@ -1,8 +1,7 @@
 function normalize_rdm(ρ::ITensor)
-    dtype = datatype(ρ)
     tr_ρ = copy(ρ)
     for i in inds(ρ; plev = 0)
-        tr_ρ *= adapt(dtype)(delta(i, prime(i)))
+        tr_ρ *= delta(scalartype(ρ), dag(i), prime(i))
     end
     return ρ / scalar(tr_ρ)
 end
diff --git a/src/sampling.jl b/src/sampling.jl
index 1051fea..635d447 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -20,7 +20,7 @@ function sample(
         bit_string = Dictionary{keytype(vertices(ψ)), Int}()
         for v in vertices(ψ)
             tensors = incoming_messages(projected_bp_cache, v)
-            ψv, ψv_dag = network(projected_bp_cache)[v], dag(prime(network(projected_bp_cache)[v]))
+            ψv, ψv_dag = network(projected_bp_cache)[v], bra_tensor(network(projected_bp_cache), v)
             push!(tensors, ψv, ψv_dag)
             seq = contraction_sequence(tensors; alg = "optimal")
             ρ = ITensors.contract(tensors; sequence = seq)
@@ -32,7 +32,7 @@ function sample(
             # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself
             set!(bit_string, v, config - 1)
             s_ind = inds(ρ)[findfirst(i -> plev(i) == 0, inds(ρ))]
-            P = adapt_like(ρ, onehot(s_ind => config))
+            P = adapt_like(ρ, dag(onehot(s_ind => config)))
             setindex_preserve!(projected_bp_cache, ψv * P, v)
 
             if v != last(vertices(ψ))
@@ -198,7 +198,13 @@ function get_one_sample(
             es = sorted_edges(norm_bmps_cache, pe)
 
             for (i, e) in enumerate(es)
-                setmessage!(norm_bmps_cache, e, ITensor[outgoing_mps[i], prime(dag(outgoing_mps[i]))])
+                mt = outgoing_mps[i]
+                # Legs unique to `mt` (not facing the network, not MPS chain bonds) are
+                # dangling non-physical legs the applied row carried into the message
+                # (e.g. charge legs of projected sites); the bra copy pairs them.
+                net_inds = virtualinds(network(norm_bmps_cache), e)
+                aux = uniqueinds(mt, cat_inds(net_inds, (inds(m) for m in outgoing_mps if m !== mt)...))
+                setmessage!(norm_bmps_cache, e, ITensor[mt, bra_tensor(mt, aux)])
             end
 
             incoming_mps = outgoing_mps
@@ -228,7 +234,7 @@ function sample_partition!(
         !isnothing(prev_v) && update_partition!(norm_bmps_cache, [NamedEdge(prev_v => v)])
         incoming_ms = incoming_messages(norm_bmps_cache, [v])
         ψv = network(norm_bmps_cache)[v]
-        ψvdag = dag(prime(ψv))
+        ψvdag = bra_tensor(network(norm_bmps_cache), v)
         ts = [incoming_ms; [ψv, ψvdag]]
         seq = contraction_sequence(ts; alg = "optimal")
         ρ = contract(ts; sequence = seq)
@@ -240,7 +246,7 @@ function sample_partition!(
         # config is 1,2,...,d, but we want 0,1...,d-1 for the sample itself
         set!(bit_string, v, config - 1)
         s_ind = inds(ρ)[findfirst(i -> plev(i) == 0, inds(ρ))]
-        P = adapt_like(ρ, onehot(s_ind => config))
+        P = adapt_like(ρ, dag(onehot(s_ind => config)))
         q = ρ_diag[config]
         logq += log(q)
         Pψv = copy(network(norm_bmps_cache)[v]) * inv(sqrt(q)) * P
@@ -268,7 +274,7 @@ function certify_sample(
     s = siteinds(ψ)
     qv = sqrt(exp(inv(oftype(logq, length(vertices(ψ)))) * logq))
     for v in vertices(ψ)
-        P = adapt_like(ψproj[v], onehot(only(s[v]) => bitstring[v] + 1))
+        P = adapt_like(ψproj[v], dag(onehot(only(s[v]) => bitstring[v] + 1)))
         setindex_preserve!(ψproj, ψproj[v] * P * inv(qv), v)
     end
 

From 0eb67bbbb4a34e63a5a2b1f3e15f6146fce33daf Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Sun, 5 Jul 2026 19:32:32 -0400
Subject: [PATCH 03/14] Support symmetric backends in loop corrections and
 re-include the example tests

The loop-correction weights now match bond indices by name (a graded tensor carries the dual of its message's axis, so full Index comparisons never match), build the antiprojector identity from the actual message legs, and record dangling aux legs before the bond relabeling so the bra construction does not misread the relabeled bond as a charge leg. The compat layer gains the legacy `itensor` reshape semantics and a shared charged-state construction for state vectors, used by `onehot` and the product-state constructor. The example scripts run on the new backend and `test_examples.jl` is back in the test suite.
---
 examples/boundarymps.jl                       |  1 +
 examples/heavyhexIsing_dynamics.jl            |  1 +
 .../hexagonal_heisenbergmodel_thermalstate.jl |  2 +-
 examples/loopcorrections.jl                   |  2 +
 src/ITensorsITensorBaseCompat/itensors.jl     | 27 ++++++------
 src/ITensorsITensorBaseCompat/ops.jl          |  7 +++-
 src/MessagePassing/loopcorrection.jl          | 41 ++++++++++++++-----
 src/TensorNetworks/tensornetworkstate.jl      |  9 ++--
 test/runtests.jl                              |  1 +
 9 files changed, 64 insertions(+), 27 deletions(-)

diff --git a/examples/boundarymps.jl b/examples/boundarymps.jl
index 4b45545..f4ba696 100644
--- a/examples/boundarymps.jl
+++ b/examples/boundarymps.jl
@@ -1,4 +1,5 @@
 using TensorNetworkQuantumSimulator
+using TensorNetworkQuantumSimulator: ITensors
 
 using Random
 Random.seed!(1634)
diff --git a/examples/heavyhexIsing_dynamics.jl b/examples/heavyhexIsing_dynamics.jl
index e6a02c7..62b87dc 100644
--- a/examples/heavyhexIsing_dynamics.jl
+++ b/examples/heavyhexIsing_dynamics.jl
@@ -1,5 +1,6 @@
 using TensorNetworkQuantumSimulator
 using Statistics
+using TensorNetworkQuantumSimulator: ITensor, ITensors
 
 function main()
     #Define the lattice
diff --git a/examples/hexagonal_heisenbergmodel_thermalstate.jl b/examples/hexagonal_heisenbergmodel_thermalstate.jl
index 3a1c76d..0c9f08e 100644
--- a/examples/hexagonal_heisenbergmodel_thermalstate.jl
+++ b/examples/hexagonal_heisenbergmodel_thermalstate.jl
@@ -1,6 +1,6 @@
 using TensorNetworkQuantumSimulator
 using TensorNetworkQuantumSimulator: scalar_factors_quotient, TensorNetworkQuantumSimulator, freenergy
-using ITensors: ITensors, ITensor
+using TensorNetworkQuantumSimulator: ITensors, ITensor
 
 function main()
     χ = 32
diff --git a/examples/loopcorrections.jl b/examples/loopcorrections.jl
index e908f0a..e5f65fd 100644
--- a/examples/loopcorrections.jl
+++ b/examples/loopcorrections.jl
@@ -1,4 +1,6 @@
 using TensorNetworkQuantumSimulator
+using TensorNetworkQuantumSimulator: ITensors
+
 using LinearAlgebra: norm
 
 using Random
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 933be68..1361b32 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -162,9 +162,17 @@ end
 # Legacy `itensor(array, inds)`: inherit the index spaces. NB: `ITensor(array, inds)`
 # with raw `Index` objects is intentionally NOT supported by ITensorBase (the space
 # is underdefined); use the indexing form, which inherits the indices' spaces, or
-# `ITensor(array, name.(inds))` to take the space from the array.
-itensor(array, is) = array[is...]
-itensor(array, is...) = array[is...]
+# `ITensor(array, name.(inds))` to take the space from the array. Like the legacy
+# constructor, a matching total length is accepted and reshaped to the index
+# dimensions (e.g. a `d^2 × d^2` two-site gate matrix over four site legs).
+function itensor(array, is...)
+    length(array) == prod(length, is) ||
+        error(
+        "array with $(length(array)) elements cannot fill indices of dimensions $(length.(is))"
+    )
+    return reshape(array, map(length, is))[is...]
+end
+itensor(array, is::Union{Tuple, AbstractVector}) = itensor(array, is...)
 
 # TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): adds a rank-0
 # `ITensor(x::Number)` constructor, which ITensorBase deliberately omits and does not plan to
@@ -254,18 +262,13 @@ function tr(t::AbstractITensor)
     return scalar(out)
 end
 
-# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), built via checked
-# `project` so it follows the index's backend. When the position's sector is charged
-# under a graded `i` the vector can't live over `i` alone, so carry the charge on a
-# derived length-1 auxiliary leg — the same construction as a charged `state` (`ops.jl`).
+# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through the
+# vector-state constructor (`state` in `ops.jl`), which follows the index's backend
+# and carries any charge on a derived auxiliary leg.
 function onehot(eltype::Type, (i, p)::Pair{<:Index})
     v = zeros(eltype, length(i))
     v[p] = one(eltype)
-    ψ = TensorAlgebra.tryproject(v, (i,))
-    isnothing(ψ) || return ψ
-    raw = TensorAlgebra.project(reshape(v, (length(v), 1)), (unnamed(i),), ())
-    aux = Index(TensorAlgebra.axes(raw, 2))
-    return nameddims(raw, (ITensorBase.name(i), ITensorBase.name(aux)))
+    return state(v, i)
 end
 onehot(p::Pair{<:Index}) = onehot(Float64, p)
 
diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/ITensorsITensorBaseCompat/ops.jl
index c5f35a0..226657a 100644
--- a/src/ITensorsITensorBaseCompat/ops.jl
+++ b/src/ITensorsITensorBaseCompat/ops.jl
@@ -30,9 +30,14 @@ function state(name::AbstractString, i::Index)
             "unknown single-site state \"$name\" (vendored states: $(sort(collect(keys(_STATE_VECTORS)))))"
         )
     end
+    return state(v, i)
+end
+# Vector form (legacy `ITensor(v, i)` for a state vector): the named state vector as
+# an `ITensor` over `i`, via checked `project` so the index axis selects the backend.
+function state(v::AbstractVector{<:Number}, i::Index)
     length(v) == length(i) ||
         error(
-        "state \"$name\" has dimension $(length(v)) but the site index has dimension $(length(i))"
+        "state vector has dimension $(length(v)) but the site index has dimension $(length(i))"
     )
     ψ = tryproject(v, (i,))
     isnothing(ψ) || return ψ
diff --git a/src/MessagePassing/loopcorrection.jl b/src/MessagePassing/loopcorrection.jl
index a37a59c..dd30010 100644
--- a/src/MessagePassing/loopcorrection.jl
+++ b/src/MessagePassing/loopcorrection.jl
@@ -19,6 +19,14 @@ end
 function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg)
     bpc = copy(bpc)
     vs = collect(vertices(eg))
+    # Aux (dangling non-physical) legs recorded before the bond relabeling below: the
+    # relabeled bond dangles in the modified network, and the live classification in
+    # `norm_factors` would misread it as a charge leg and unprime it on the bra.
+    aux = if network(bpc) isa TensorNetworkState
+        Dictionary(vs, [auxinds(network(bpc), v) for v in vs])
+    else
+        nothing
+    end
     es =
         unique(collect(Iterators.flatten(boundary_edges(bpc, [v]; dir = :out) for v in vs)))
     updated_es = NamedEdge[]
@@ -35,28 +43,37 @@ function sim_edgeinduced_subgraph(bpc::BeliefPropagationCache, eg)
             ms = messages(bpc)
             set!(ms, reverse(e), mer)
             t = network(bpc)[src(e)]
-            t_inds = filter(i -> i ∈ linds, inds(t))
+            # Match by name: on a graded backend the network tensor carries the dual of
+            # the message's axis, so a full `Index` comparison never matches.
+            t_inds = commoninds(t, linds)
             if !isempty(t_inds)
                 t_ind = only(t_inds)
-                t_ind_pos = findfirst(x -> x == t_ind, linds)
+                t_ind_pos = findfirst(x -> name(x) == name(t_ind), linds)
                 t = replaceind(t, t_ind, linds_sim[t_ind_pos])
                 setindex_preserve!(bpc, t, src(e))
             end
             push!(updated_es, e)
 
             if e ∈ edges(eg) || reverse(e) ∈ edges(eg)
-                row_inds, col_inds = linds, linds_sim
+                me = message(bpc, e)
+                # The identity legs take their axes from the actual message legs (`me`
+                # for the rows, the relabeled `mer` for the columns) so `ap - me * mer`
+                # lines up on every backend: the two messages carry mutually dual copies
+                # of the bond axes. The domain of the fused identity comes out dualized
+                # relative to the passed indices, so the columns go in `dag`ed.
+                row_inds = Index[only(commoninds(me, [l])) for l in linds]
+                col_inds = Index[only(commoninds(mer, [l])) for l in linds_sim]
                 if network(bpc) isa TensorNetworkState
-                    row_inds = vcat(row_inds, dag.(prime.(row_inds)))
-                    col_inds = vcat(col_inds, dag.(prime.(col_inds)))
+                    append!(row_inds, Index[only(commoninds(me, [prime(l)])) for l in linds])
+                    append!(col_inds, Index[only(commoninds(mer, [prime(l)])) for l in linds_sim])
                 end
-                ap = adapt_like(message(bpc, e), identity_tensor(row_inds, col_inds))
-                ap = ap - message(bpc, e) * mer
+                ap = adapt_like(me, identity_tensor(row_inds, dag.(col_inds)))
+                ap = ap - me * mer
                 push!(antiprojectors, ap)
             end
         end
     end
-    return bpc, antiprojectors
+    return bpc, antiprojectors, aux
 end
 
 #Get the all edges incident to the region specified by the vector of edges passed
@@ -78,10 +95,14 @@ end
 function weight(bpc::BeliefPropagationCache, eg)
     vs = collect(vertices(eg))
     es = collect(edges(eg))
-    bpc, antiprojectors = sim_edgeinduced_subgraph(bpc, eg)
+    bpc, antiprojectors, aux = sim_edgeinduced_subgraph(bpc, eg)
     incoming_ms =
         ITensor[message(bpc, e) for e in boundary_edges(bpc, es)]
-    local_tensors = collect(Iterators.flatten(bp_factors(bpc, v) for v in vs))
+    local_tensors = if isnothing(aux)
+        collect(Iterators.flatten(bp_factors(bpc, v) for v in vs))
+    else
+        collect(Iterators.flatten(norm_factors(network(bpc), [v]; auxinds_f = u -> aux[u]) for v in vs))
+    end
     ts = [incoming_ms; local_tensors; antiprojectors]
     seq = any(hasqns.(ts)) ? contraction_sequence(ts; alg = "optimal") : contraction_sequence(ts; alg = "omeinsum", optimizer = GreedyMethod())
     return scalar(contract(ts; sequence = seq))
diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl
index 75cbf83..b4140ed 100644
--- a/src/TensorNetworks/tensornetworkstate.jl
+++ b/src/TensorNetworks/tensornetworkstate.jl
@@ -67,12 +67,15 @@ bra_tensor(tns::TensorNetworkState, v) = bra_tensor(tns[v], auxinds(tns, v))
 # The dangling non-physical legs of a vertex tensor: dangling legs that are not site indices.
 auxinds(tns::TensorNetworkState, v) = Index[i for i in setdiff(uniqueinds(tns, v), siteinds(tns, v))]
 
-function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I")
+# `auxinds_f` overrides the live aux-leg classification. The loop-correction weights
+# need this: they deliberately relabel a bond so it dangles in the modified network,
+# and the live classification would misread it as a charge leg.
+function norm_factors(tns::TensorNetworkState, verts::Vector; op_strings::Function = v -> "I", auxinds_f::Function = v -> auxinds(tns, v))
     factors = ITensor[]
     for v in verts
         sinds = siteinds(tns, v)
         tnv = tns[v]
-        tnv_dag = bra_tensor(tns, v)
+        tnv_dag = bra_tensor(tnv, auxinds_f(v))
         if op_strings(v) == "ρ" || isempty(sinds)
             append!(factors, ITensor[tnv, tnv_dag])
         elseif op_strings(v) == "I"
@@ -170,7 +173,7 @@ function tensornetworkstate(eltype, f::Function, g::AbstractGraph, siteinds::Dic
         if tnv isa String
             set!(tensors, v, adapt(eltype)(ITensors.state(f(v), only(siteinds[v]))))
         elseif tnv isa Vector{<:Number}
-            set!(tensors, v, adapt(eltype)(ITensors.ITensor(f(v), only(siteinds[v]))))
+            set!(tensors, v, adapt(eltype)(ITensors.state(f(v), only(siteinds[v]))))
         else
             error("Unrecognized local state constructor. Currently supported: Strings and Vectors.")
         end
diff --git a/test/runtests.jl b/test/runtests.jl
index 641e25a..1cec2bc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,4 +11,5 @@ using Test
     include("test_sampling.jl")
     include("test_truncate.jl")
     include("test_contraction_sequences.jl")
+    include("test_examples.jl")
 end

From 21d5aca82106c9cf1e6a76e70a82d4f6165a8a74 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Sun, 5 Jul 2026 21:31:00 -0400
Subject: [PATCH 04/14] Rewrite symmetric_gauge in the BP perspective with
 rename-based absorption

The gauge transformation now works directly on the BP cache: per edge, take Hermitian square roots and pseudo-inverse roots of the two messages with sqrth_invsqrth_safe, absorb them into the endpoint tensors by renaming the state's bond indices onto the roots' domain names, and SVD the bond matrix formed by contracting the two roots' codomain legs. The factorization directions and names line up so no conjugation or arrow flips are needed on any backend (dense, graded, TensorKit). Small eigenvalues are clamped to zero through the pseudo-inverse instead of the legacy regularization shift.

The safe spectral functions require Hermitian input, so this and the other square-root callsites (simple_update, entanglement) project with project_hermitian first. The compat sim now rebuilds the index name keeping tags and prime level, since uniquename mints bare names.

Co-authored-by: Claude <noreply@anthropic.com>
---
 src/Apply/simple_update.jl                |  5 +-
 src/ITensorsITensorBaseCompat/itensors.jl |  9 ++-
 src/entanglement.jl                       |  3 +-
 src/imports.jl                            |  1 +
 src/symmetric_gauge.jl                    | 82 +++++++++++++----------
 5 files changed, 59 insertions(+), 41 deletions(-)

diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl
index d7aa292..3fc5760 100644
--- a/src/Apply/simple_update.jl
+++ b/src/Apply/simple_update.jl
@@ -35,8 +35,11 @@ function simple_update(
         envs_v2 = filter(env -> hascommoninds(env, ψ⃗[2]), envs)
         @assert all(ndims(env) == 2 for env in vcat(envs_v1, envs_v2))
 
+        # The environments are hermitian only up to numerical noise, so project before
+        # the square roots (which require hermitian input).
         sqrt_invsqrt = env -> sqrth_invsqrth_safe(
-            env, (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0
+            project_hermitian(env, (inds(env)[1],), (inds(env)[2],)),
+            (inds(env)[1],), (inds(env)[2],); atol = sqrt_cutoff, rtol = 0
         )
         sqrt_inv_sqrt_envs_v1 = map(sqrt_invsqrt, envs_v1)
         sqrt_inv_sqrt_envs_v2 = map(sqrt_invsqrt, envs_v2)
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 1361b32..dd50127 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -91,8 +91,13 @@ dim(is::Union{Tuple, AbstractVector}) = isempty(is) ? 1 : prod(length, is)
 #
 # Index operations.
 #
-# A fresh index with the same length (legacy `sim`).
-sim(i::Index) = ITensorBase.uniquename(i)
+# A fresh index with the same length, tags, and prime level (legacy `sim`).
+# `uniquename` mints bare names, so rebuild the name keeping the decoration.
+function sim(i::Index)
+    n = ITensorBase.name(i)
+    n_sim = ITensorBase.IndexName(; tags = ITensorBase.tags(n), plev = ITensorBase.plev(n))
+    return ITensorBase.named(ITensorBase.unnamed(i), n_sim)
+end
 sim(is::Union{Tuple, AbstractVector{<:Index}}) = map(sim, is)
 
 # Conjugate (legacy `dag`): `conj` the tensor, and on bare indices flip the sector arrows
diff --git a/src/entanglement.jl b/src/entanglement.jl
index c766271..25496da 100644
--- a/src/entanglement.jl
+++ b/src/entanglement.jl
@@ -81,7 +81,8 @@ function renyi_entropy(
     m1, m2 = message(bp_cache, e), message(bp_cache, reverse(e))
     edge_ind = only(virtualinds(bp_cache, e))
     root_m2 = sqrth_safe(
-        m2, (inds(m2)[1],), (inds(m2)[2],);
+        project_hermitian(m2, (inds(m2)[1],), (inds(m2)[2],)),
+        (inds(m2)[1],), (inds(m2)[2],);
         atol = 10 * eps(real(scalartype(m2))), rtol = 0
     )
 
diff --git a/src/imports.jl b/src/imports.jl
index 5a58c24..6189c83 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -69,6 +69,7 @@ using .ITensorsITensorBaseCompat:
 using ITensorBase: ITensorBase, Index, ITensor, name, plev, tags, unnamed
 using TensorAlgebra: trivialrange
 using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe
+using MatrixAlgebraKit: project_hermitian
 
 using Adapt: adapt
 
diff --git a/src/symmetric_gauge.jl b/src/symmetric_gauge.jl
index f5aa129..7c2bda5 100644
--- a/src/symmetric_gauge.jl
+++ b/src/symmetric_gauge.jl
@@ -1,5 +1,4 @@
 function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10 * eps(real(scalartype(bp_cache))), kwargs...)
-    dtype = datatype(bp_cache)
     tn = network(bp_cache)
     !(tn isa TensorNetworkState) && error("Can only transform TensorNetworkStates to the symmetric gauge")
     for e in edges(tn)
@@ -7,49 +6,58 @@ function symmetric_gauge!(bp_cache::BeliefPropagationCache; regularization = 10
         ψvsrc, ψvdst = tn[vsrc], tn[vdst]
 
         edge_ind = commoninds(ψvsrc, ψvdst)
+        edge_ind_p = prime.(edge_ind)
         edge_ind_sim = sim(edge_ind)
 
-        X_D, X_U = safe_eigen(message(bp_cache, e); ishermitian = true, cutoff = nothing)
-        Y_D, Y_U = safe_eigen(message(bp_cache, reverse(e)); ishermitian = true, cutoff = nothing)
-        X_D, Y_D = ITensors.map_diag(x -> x + regularization, X_D),
-            ITensors.map_diag(x -> x + regularization, Y_D)
-
-        rootX_D, rootY_D = ITensors.map_diag(x -> sqrt(x), X_D), ITensors.map_diag(x -> sqrt(x), Y_D)
-        inv_rootX_D, inv_rootY_D = ITensors.map_diag(x -> inv(sqrt(x)), X_D),
-            ITensors.map_diag(x -> inv(sqrt(x)), Y_D)
-        rootX = X_U * rootX_D * prime(dag(X_U))
-        rootY = Y_U * rootY_D * prime(dag(Y_U))
-        inv_rootX = X_U * inv_rootX_D * prime(dag(X_U))
-        inv_rootY = Y_U * inv_rootY_D * prime(dag(Y_U))
-
-        ψvsrc, ψvdst = noprime(ψvsrc * inv_rootX), noprime(ψvdst * inv_rootY)
-
-        Ce = rootX
-        Ce = Ce * replaceinds(rootY, edge_ind, edge_ind_sim)
-
-        U, S, V = svd(Ce, edge_ind; kwargs...)
-
-        new_edge_ind = Index[Index(dim(commoninds(S, U)), tags(first(edge_ind)))]
-
-        ψvsrc = replaceinds(ψvsrc * U, commoninds(S, U), new_edge_ind)
-        ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_sim)
-        ψvdst = replaceinds(ψvdst * V, commoninds(V, S), new_edge_ind)
-
-
-        S = replaceinds(
-            S,
-            [commoninds(S, U)..., commoninds(S, V)...] =>
-                [new_edge_ind..., prime(new_edge_ind)...],
+        # Hermitian square roots (and pseudo-inverse roots) of the two messages, as
+        # operators with codomain on the unprimed and domain on the primed copies of the
+        # bond indices. The messages are hermitian only up to numerical noise, so project
+        # first. Eigenvalues below `regularization` are clamped to zero.
+        me = project_hermitian(message(bp_cache, e), Tuple(edge_ind), Tuple(edge_ind_p))
+        mer = project_hermitian(message(bp_cache, reverse(e)), Tuple(edge_ind), Tuple(edge_ind_p))
+        rootX, inv_rootX = sqrth_invsqrth_safe(
+            me, Tuple(edge_ind), Tuple(edge_ind_p);
+            atol = regularization, rtol = 0,
+        )
+        rootY, inv_rootY = sqrth_invsqrth_safe(
+            mer, Tuple(edge_ind), Tuple(edge_ind_p);
+            atol = regularization, rtol = 0,
         )
 
-        sqrtS = ITensors.map_diag(sqrt, S)
-        ψvsrc = noprime(ψvsrc * sqrtS)
-        ψvdst = noprime(ψvdst * sqrtS)
+        # Absorb the inverse roots through their domain legs: move the state's bond
+        # indices onto the primed names so each pairs with its dual copy, and the
+        # result comes back out on the unprimed (codomain) names.
+        ψvsrc = replaceinds(ψvsrc, edge_ind, edge_ind_p) * inv_rootX
+        ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_p) * inv_rootY
+
+        # Bond matrix: contract the codomain legs of the two roots (one from each side
+        # of the edge, mutually dual). Its open legs are then dual to the state
+        # tensors' bond legs, so the SVD factors below absorb without any flips.
+        Ce = rootX * replaceinds(rootY, edge_ind_p, edge_ind_sim)
+
+        U, S, V = svd(Ce, edge_ind_p; kwargs...)
+        u, v = commoninds(S, U), commoninds(S, V)
+
+        ψvsrc = replaceinds(ψvsrc, edge_ind, edge_ind_p) * U
+        ψvdst = replaceinds(ψvdst, edge_ind, edge_ind_sim) * V
+
+        # Split the singular values symmetrically into both endpoints. `S`'s legs are
+        # dual to the tensors' new bond legs on both sides, so these contract directly.
+        # `S` is diagonal, so this hits the diagonal fast path on every backend.
+        sqrtS = sqrth_safe(S, Tuple(u), Tuple(v); atol = 0, rtol = 0)
+        ψvsrc = ψvsrc * sqrtS
+        ψvdst = ψvdst * sqrtS
+
+        new_edge_ind = Index[ITensors.settags(only(v), tags(first(edge_ind)))]
+        ψvsrc = replaceinds(ψvsrc, v, new_edge_ind)
+        ψvdst = replaceinds(ψvdst, u, new_edge_ind)
         setindex_preserve!(bp_cache, ψvsrc, vsrc)
         setindex_preserve!(bp_cache, ψvdst, vdst)
 
-        setmessage!(bp_cache, e, S)
-        setmessage!(bp_cache, reverse(e), dag(S))
+        # The gauged network's messages are the singular values, with the unprimed leg
+        # carrying the producing side's bond copy (the message convention).
+        setmessage!(bp_cache, e, replaceinds(S, [only(u), only(v)], [prime(only(new_edge_ind)), only(new_edge_ind)]))
+        setmessage!(bp_cache, reverse(e), replaceinds(S, [only(u), only(v)], [only(new_edge_ind), prime(only(new_edge_ind))]))
     end
 
     return bp_cache

From 406a5ab24772ba8eed2dfc3c059bf30d3fecd04d Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 6 Jul 2026 11:14:11 -0400
Subject: [PATCH 05/14] Pin the ITensorBase and TensorAlgebra branches while
 the symmetry support is in review

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 1aa3886..9f01eb3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -54,5 +54,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 test = ["Test", "Random"]
 
 [sources]
-ITensorBase = {path = "/Users/mfishman/.julia/dev/ITensorBase/.worktrees/mf/tensoralgebra-0.17"}
-TensorAlgebra = {path = "/Users/mfishman/.julia/dev/TensorAlgebra/.worktrees/mf/project-trailing-axes"}
+ITensorBase = {url = "https://github.com/ITensor/ITensorBase.jl", rev = "mf/tensoralgebra-0.17"}
+TensorAlgebra = {url = "https://github.com/ITensor/TensorAlgebra.jl", rev = "mf/project-trailing-axes"}

From c8eba66a582cf002805b0231e91a4152293114c5 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Mon, 6 Jul 2026 13:12:43 -0400
Subject: [PATCH 06/14] Tidy the ITensors compat: sim, onehot, and the identity
 map

`sim` is now `uniquename` on the index, which keeps the tags and prime level. `onehot` and the state-vector constructor share a `project_aux` helper that carries a charge on a derived auxiliary leg only when the basis vector needs it. The order-2 identity map is spelled `id(codomain, domain)`, following TensorKit, and the direct identity-map callers (`tr`, `normalize_rdm`, `bilinearform`) use it, while `delta` keeps its own order-2 method for the count-dispatched copy-tensor cases. Factorization and message bond names mint bare `IndexName`s via the name type.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/Forms/bilinearform.jl                     |  2 +-
 .../ITensorsITensorBaseCompat.jl              |  5 +-
 src/ITensorsITensorBaseCompat/itensors.jl     | 54 ++++++++++---------
 src/ITensorsITensorBaseCompat/ops.jl          | 20 ++++---
 src/MessagePassing/beliefpropagationcache.jl  |  4 +-
 src/entanglement.jl                           |  4 +-
 src/imports.jl                                |  2 +-
 src/rdm.jl                                    |  2 +-
 8 files changed, 51 insertions(+), 42 deletions(-)

diff --git a/src/Forms/bilinearform.jl b/src/Forms/bilinearform.jl
index 7c397d9..234f6aa 100644
--- a/src/Forms/bilinearform.jl
+++ b/src/Forms/bilinearform.jl
@@ -18,7 +18,7 @@ function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState)
     sinds = siteinds(ket)
     verts = collect(vertices(ket))
     bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts]))
-    operator_tensors = [reduce(*, ITensor[delta(scalartype(ket), dag(sind), prime(sind)) for sind in sinds[v]]) for v in verts]
+    operator_tensors = [reduce(*, ITensor[id(scalartype(ket), (dag(sind),), (dag(prime(sind)),)) for sind in sinds[v]]) for v in verts]
     operator = TensorNetworkState(Dictionary(verts, operator_tensors))
     return BilinearForm(ket, operator, bra)
 end
diff --git a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
index b32f293..72a0210 100644
--- a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
+++ b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
@@ -26,11 +26,12 @@ include("ops.jl")
 
 export
     # Index access and set algebra
-    inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds,
+    inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds,
+    hascommoninds,
     # Index operations
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
     # ITensor construction
-    itensor, random_itensor, scalar, delta, onehot,
+    itensor, random_itensor, scalar, delta, id, onehot,
     # Factorizations
     qr, svd, eigen, factorize, factorize_svd,
     # Diagonal manipulation
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index dd50127..b9a4fcd 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -92,12 +92,7 @@ dim(is::Union{Tuple, AbstractVector}) = isempty(is) ? 1 : prod(length, is)
 # Index operations.
 #
 # A fresh index with the same length, tags, and prime level (legacy `sim`).
-# `uniquename` mints bare names, so rebuild the name keeping the decoration.
-function sim(i::Index)
-    n = ITensorBase.name(i)
-    n_sim = ITensorBase.IndexName(; tags = ITensorBase.tags(n), plev = ITensorBase.plev(n))
-    return ITensorBase.named(ITensorBase.unnamed(i), n_sim)
-end
+sim(i::Index) = ITensorBase.uniquename(i)
 sim(is::Union{Tuple, AbstractVector{<:Index}}) = map(sim, is)
 
 # Conjugate (legacy `dag`): `conj` the tensor, and on bare indices flip the sector arrows
@@ -230,20 +225,29 @@ function diagonaltensor(
     return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is))
 end
 
+# The identity map between the codomain and domain index groups (TensorKit's `id`), built
+# via checked `project` on a dense identity so the index axes select the backend (dense,
+# graded, `TensorMap`). The result carries exactly the given indices with their own axes
+# (`project` dualizes the domain axes). Callers arrow the indices, as for `identity_tensor`:
+# the fused codomain and domain must be arrow-opposite (e.g. `id(elt, (i,), (dag(j),))`); an
+# identity between incompatible gradings throws.
+function id(eltype::Type, codomain::Tuple, domain::Tuple)
+    m = Matrix{eltype}(LinearAlgebra.I, prod(length, codomain), prod(length, domain))
+    return TensorAlgebra.project(m, codomain, domain)
+end
+
 delta(eltype::Type, is::Tuple) = diagonaltensor(ones(eltype, minimum(length, is)), is)
-# The order-2 delta is the identity map over the index pair, built via checked
-# `project` so the index axes select the backend (dense, graded, `TensorMap`). The
-# result carries exactly the given indices with their own axes (`project` dualizes
-# domain axes, so `dag(j)` in the domain slot passes `j`'s axis through) — the legacy
-# QN-`delta` convention, where callers arrow the indices themselves. The pair must be
-# arrow-opposite (axis of `j` the dual of axis of `i`, e.g. `(dag(i), prime(i))`); an
-# identity between incompatible gradings throws. Higher-order deltas stay on the dense
-# `diagonaltensor` path above (a super-diagonal generally cannot be embedded while
-# preserving a nontrivial symmetry).
+# The order-2 delta is the identity map over the index pair, built via checked `project` so
+# the index axes select the backend (dense, graded, `TensorMap`). Legacy QN-`delta`
+# convention: the pair is arrow-opposite and the caller passes `(i, j)` with `j` un-dualized.
+# It coincides with `id` at order 2, but stays its own method (its own calling convention,
+# and the count dispatch below routes graded default messages here). Higher-order deltas stay
+# on the dense `diagonaltensor` path above (a super-diagonal generally cannot be embedded
+# while preserving a nontrivial symmetry).
 function delta(eltype::Type, is::Tuple{Index, Index})
     i, j = is
-    id = Matrix{eltype}(LinearAlgebra.I, length(i), length(j))
-    return TensorAlgebra.project(id, (i,), (dag(j),))
+    m = Matrix{eltype}(LinearAlgebra.I, length(i), length(j))
+    return TensorAlgebra.project(m, (i,), (dag(j),))
 end
 delta(eltype::Type, is::Index...) = delta(eltype, is)
 delta(eltype::Type, is::AbstractVector{<:Index}) = delta(eltype, Tuple(is))
@@ -260,20 +264,20 @@ delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is))
 function tr(t::AbstractITensor)
     out = copy(t)
     for i in inds(t; plev = 0)
-        # The delta legs carry the duals of `t`'s prime pair (`i` and its primed dual),
-        # so each leg contracts its partner on any backend.
-        out *= delta(scalartype(t), dag(i), prime(i))
+        # The identity-map legs carry the duals of `t`'s prime pair (`i` and its primed
+        # dual), so each leg contracts its partner on any backend.
+        out *= id(scalartype(t), (dag(i),), (dag(prime(i)),))
     end
     return scalar(out)
 end
 
-# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through the
-# vector-state constructor (`state` in `ops.jl`), which follows the index's backend
-# and carries any charge on a derived auxiliary leg.
+# One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through `project_aux`
+# (in `ops.jl`), which follows the index's backend and carries the basis vector's charge on
+# a derived auxiliary leg.
 function onehot(eltype::Type, (i, p)::Pair{<:Index})
     v = zeros(eltype, length(i))
     v[p] = one(eltype)
-    return state(v, i)
+    return project_aux(v, i)
 end
 onehot(p::Pair{<:Index}) = onehot(Float64, p)
 
@@ -509,7 +513,7 @@ end
 #
 # Index fusion. The legacy `combiner` is retired (no compat shim); call sites fuse
 # index groups with the next-gen `matricize(t, row_inds => row_name, col_inds =>
-# col_name)` (minting each fused name via `name(uniquename(i))`) or build a fused
+# col_name)` (minting each fused name via `uniquename(IndexName)`) or build a fused
 # identity with `Base.one`, both of which are graded-capable. `matricize` is the
 # `TensorAlgebra` generic, extended by ITensorBase for named tensors.
 using TensorAlgebra: matricize
diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/ITensorsITensorBaseCompat/ops.jl
index 226657a..14ee995 100644
--- a/src/ITensorsITensorBaseCompat/ops.jl
+++ b/src/ITensorsITensorBaseCompat/ops.jl
@@ -32,24 +32,28 @@ function state(name::AbstractString, i::Index)
     end
     return state(v, i)
 end
-# Vector form (legacy `ITensor(v, i)` for a state vector): the named state vector as
-# an `ITensor` over `i`, via checked `project` so the index axis selects the backend.
-function state(v::AbstractVector{<:Number}, i::Index)
+# Project a raw vector as a state `ITensor` over `i`, adding an auxiliary leg only when the
+# vector cannot live in the flux-zero space over `i` alone. The index axis selects the
+# backend (dense, graded, `TensorMap`). Shared by the state constructors and `onehot`.
+function project_aux(v::AbstractVector{<:Number}, i::Index)
     length(v) == length(i) ||
         error(
         "state vector has dimension $(length(v)) but the site index has dimension $(length(i))"
     )
     ψ = tryproject(v, (i,))
     isnothing(ψ) || return ψ
-    # The state carries a charge under the site index's grading (e.g. "Dn" on a
-    # U(1) site), so it can't live in the flux-zero space over `i` alone. Carry the
-    # charge on an explicit length-1 auxiliary leg: project with a trailing axis so
-    # the backend derives the leg's sector from the state vector, then wrap the
-    # derived axis in a freshly named `Index`.
+    # The vector carries a charge under `i`'s grading (e.g. "Dn" on a U(1) site, or a
+    # one-hot on a graded link), so it can't live in the flux-zero space over `i` alone.
+    # Carry the charge on an explicit length-1 auxiliary leg: project with a trailing axis
+    # so the backend derives the leg's sector from the vector, then wrap the derived axis
+    # in a freshly named `Index`.
     raw = project(reshape(v, (length(v), 1)), (unnamed(i),), ())
     aux = Index(TensorAlgebra.axes(raw, 2))
     return nameddims(raw, (ITensorBase.name(i), ITensorBase.name(aux)))
 end
+# Vector form (legacy `ITensor(v, i)` for a state vector): the state vector as an `ITensor`
+# over `i`.
+state(v::AbstractVector{<:Number}, i::Index) = project_aux(v, i)
 
 #
 # Operators. Legacy ITensors exposes operators through the `OpName` / `SiteType`
diff --git a/src/MessagePassing/beliefpropagationcache.jl b/src/MessagePassing/beliefpropagationcache.jl
index 93fbc70..224d993 100644
--- a/src/MessagePassing/beliefpropagationcache.jl
+++ b/src/MessagePassing/beliefpropagationcache.jl
@@ -170,8 +170,8 @@ function loop_correlation(bpc::BeliefPropagationCache, loop::Vector{<:NamedEdge}
     seq = contraction_sequence(tensors; alg = "omeinsum", optimizer = GreedyMethod())
     t = contract(tensors; sequence = seq)
 
-    row_name = name(ITensorBase.uniquename(first(e_virtualinds)))
-    col_name = name(ITensorBase.uniquename(first(e_virtualinds_sim)))
+    row_name = ITensorBase.uniquename(ITensorBase.IndexName)
+    col_name = ITensorBase.uniquename(ITensorBase.IndexName)
     t = ITensors.matricize(t, Tuple(e_virtualinds) => row_name, Tuple(e_virtualinds_sim) => col_name)
     t = adapt(Vector{ComplexF64})(t)
     t = ITensors.array(t)
diff --git a/src/entanglement.jl b/src/entanglement.jl
index 25496da..0a93c0c 100644
--- a/src/entanglement.jl
+++ b/src/entanglement.jl
@@ -30,8 +30,8 @@ end
 
 function matricize(a::ITensor, row_inds = filter(i -> plev(i) ==0, inds(a)))
     col_inds = prime.(row_inds)
-    row_name = name(ITensorBase.uniquename(first(row_inds)))
-    col_name = name(ITensorBase.uniquename(first(col_inds)))
+    row_name = ITensorBase.uniquename(ITensorBase.IndexName)
+    col_name = ITensorBase.uniquename(ITensorBase.IndexName)
     return ITensors.array(ITensors.matricize(a, Tuple(row_inds) => row_name, Tuple(col_inds) => col_name))
 end
 
diff --git a/src/imports.jl b/src/imports.jl
index 6189c83..04ba865 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -56,7 +56,7 @@ import .ITensorsITensorBaseCompat as ITensors
 using .ITensorsITensorBaseCompat:
     inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, cat_inds,
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
-    itensor, random_itensor, scalar, delta, onehot,
+    itensor, random_itensor, scalar, delta, id, onehot,
     qr, svd, eigen, factorize, factorize_svd,
     map_diag, map_diag!,
     scalartype, datatype, array, data,
diff --git a/src/rdm.jl b/src/rdm.jl
index f07c6c4..8a40045 100644
--- a/src/rdm.jl
+++ b/src/rdm.jl
@@ -1,7 +1,7 @@
 function normalize_rdm(ρ::ITensor)
     tr_ρ = copy(ρ)
     for i in inds(ρ; plev = 0)
-        tr_ρ *= delta(scalartype(ρ), dag(i), prime(i))
+        tr_ρ *= id(scalartype(ρ), (dag(i),), (dag(prime(i)),))
     end
     return ρ / scalar(tr_ρ)
 end

From 40dca82855ee5db2ced7f1567b7a15c42503c2fc Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 11:44:33 -0400
Subject: [PATCH 07/14] Make the order-2 compat identities device-following

The order-2 identity constructions (`tr`, `normalize_rdm`, `bilinearform`) now build their identity maps through ITensorBase's device-following `Base.one` instead of a CPU-materialized dense identity routed through `project`, so the result follows the input's device (fixing a GPU regression) and stays correct on graded backends.

`tr` and `normalize_rdm` build the identity in a single `one(t, codomain, prime.(codomain))` over all prime pairs at once. The domain is `prime.(codomain)` rather than `inds(t; plev=1)`: `plev` filtering does not preserve the pairing order between the unprimed and primed groups, so pairing by prime keeps the trace independent of storage order.

`bilinearform`'s ket tensor carries each site index but not its prime, so the tensor itself can't serve as the identity's shape prototype. A new `similar_map` helper allocates an identity-shaped tensor over the given codomain and domain following the prototype's backend and element type, which `one` then fills.

Also drops the now-unused `id` (the direct identity-map callers were its only users) and the `map_diag`/`map_diag!` diagonal-map shims (no remaining callsites, and the spectral matrix functions cover those uses).
---
 src/Forms/bilinearform.jl                     |  2 +-
 .../ITensorsITensorBaseCompat.jl              |  4 +-
 src/ITensorsITensorBaseCompat/itensors.jl     | 89 ++++++-------------
 src/imports.jl                                |  3 +-
 src/rdm.jl                                    |  6 +-
 5 files changed, 32 insertions(+), 72 deletions(-)

diff --git a/src/Forms/bilinearform.jl b/src/Forms/bilinearform.jl
index 234f6aa..aa199d5 100644
--- a/src/Forms/bilinearform.jl
+++ b/src/Forms/bilinearform.jl
@@ -18,7 +18,7 @@ function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState)
     sinds = siteinds(ket)
     verts = collect(vertices(ket))
     bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts]))
-    operator_tensors = [reduce(*, ITensor[id(scalartype(ket), (dag(sind),), (dag(prime(sind)),)) for sind in sinds[v]]) for v in verts]
+    operator_tensors = [one(similar_map(ket[v], sinds[v], prime.(sinds[v])), sinds[v], prime.(sinds[v])) for v in verts]
     operator = TensorNetworkState(Dictionary(verts, operator_tensors))
     return BilinearForm(ket, operator, bra)
 end
diff --git a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
index 72a0210..21b475c 100644
--- a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
+++ b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
@@ -31,11 +31,9 @@ export
     # Index operations
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
     # ITensor construction
-    itensor, random_itensor, scalar, delta, id, onehot,
+    itensor, random_itensor, scalar, delta, onehot,
     # Factorizations
     qr, svd, eigen, factorize, factorize_svd,
-    # Diagonal manipulation
-    map_diag, map_diag!,
     # Storage / element-type accessors
     scalartype, datatype, array, data,
     # Dense / quantum-number no-ops
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index b9a4fcd..4d4d0a3 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -3,7 +3,7 @@
 #
 # Strategy (see ITensorDevelopmentPlans api_migration_map.md):
 #   - Names below are thin wrappers over ITensorBase / TensorAlgebra / MatrixAlgebraKit.
-#   - The factorization return shapes, `map_diag`, the operator/SiteType system, and
+#   - The factorization return shapes, the operator/SiteType system, and
 #     the boundary-MPS (ITensorMPS) paths are NOT wrapped here; they need callsite
 #     translation or upstream stack work and are tracked separately. The legacy
 #     `combiner` is retired rather than wrapped: call sites use the next-gen fusion
@@ -225,23 +225,28 @@ function diagonaltensor(
     return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is))
 end
 
-# The identity map between the codomain and domain index groups (TensorKit's `id`), built
-# via checked `project` on a dense identity so the index axes select the backend (dense,
-# graded, `TensorMap`). The result carries exactly the given indices with their own axes
-# (`project` dualizes the domain axes). Callers arrow the indices, as for `identity_tensor`:
-# the fused codomain and domain must be arrow-opposite (e.g. `id(elt, (i,), (dag(j),))`); an
-# identity between incompatible gradings throws.
-function id(eltype::Type, codomain::Tuple, domain::Tuple)
-    m = Matrix{eltype}(LinearAlgebra.I, prod(length, codomain), prod(length, domain))
-    return TensorAlgebra.project(m, codomain, domain)
+# Allocate an identity-map-shaped tensor with undefined data over the `codomain`/`domain`
+# indices, following `prototype`'s backend and element type with the domain axes dualized in
+# storage (the `similar_map` map convention). `prototype` only donates the backend, so unlike
+# `Base.one(a, codomain, domain)` (which needs `a` to already carry the map's legs) a caller
+# whose prototype lacks a leg still builds a device-following identity via
+# `one(similar_map(prototype, codomain, domain), codomain, domain)`.
+function similar_map(prototype::AbstractITensor, eltype::Type, codomain, domain)
+    raw = TensorAlgebra.similar_map(
+        unnamed(prototype), eltype, unnamed.(codomain), unnamed.(domain)
+    )
+    return nameddims(raw, (name.(codomain)..., name.(domain)...))
+end
+function similar_map(prototype::AbstractITensor, codomain, domain)
+    return similar_map(prototype, scalartype(prototype), codomain, domain)
 end
 
 delta(eltype::Type, is::Tuple) = diagonaltensor(ones(eltype, minimum(length, is)), is)
 # The order-2 delta is the identity map over the index pair, built via checked `project` so
 # the index axes select the backend (dense, graded, `TensorMap`). Legacy QN-`delta`
 # convention: the pair is arrow-opposite and the caller passes `(i, j)` with `j` un-dualized.
-# It coincides with `id` at order 2, but stays its own method (its own calling convention,
-# and the count dispatch below routes graded default messages here). Higher-order deltas stay
+# It coincides with the order-2 identity map, but stays its own method (its own calling
+# convention, and the count dispatch below routes graded default messages here). Higher-order deltas stay
 # on the dense `diagonaltensor` path above (a super-diagonal generally cannot be embedded
 # while preserving a nontrivial symmetry).
 function delta(eltype::Type, is::Tuple{Index, Index})
@@ -255,20 +260,17 @@ delta(is::Tuple) = delta(Float64, is)
 delta(is::Index...) = delta(Float64, is)
 delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is))
 
-# Trace of an ITensor over its prime pairs (legacy `tr`): contract each unprimed index
-# with its prime via a `delta`, the same construction `normalize_rdm` uses. This is the
-# index-paired definition — independent of storage order — rather than `sum` of the
-# underlying array's storage diagonal, which is only the trace for a rank-2 tensor whose
-# legs happen to be ordered to align. Accessed qualified (`ITensors.tr`) so it doesn't
-# shadow `LinearAlgebra.tr`, which TNQS still calls on plain matrices.
+# Trace of an ITensor over its prime pairs (legacy `tr`): contract with the identity map that
+# pairs every unprimed index with its prime, the same construction `normalize_rdm` uses. This
+# is the index-paired definition — independent of storage order — rather than `sum` of the
+# underlying array's storage diagonal, which is only the trace for a rank-2 tensor whose legs
+# happen to be ordered to align. The domain is `prime.(codomain)` (not `inds(t; plev=1)`) so
+# each codomain index pairs with its own prime; `plev` filtering does not preserve that
+# pairing order. Accessed qualified (`ITensors.tr`) so it doesn't shadow `LinearAlgebra.tr`,
+# which TNQS still calls on plain matrices.
 function tr(t::AbstractITensor)
-    out = copy(t)
-    for i in inds(t; plev = 0)
-        # The identity-map legs carry the duals of `t`'s prime pair (`i` and its primed
-        # dual), so each leg contracts its partner on any backend.
-        out *= id(scalartype(t), (dag(i),), (dag(prime(i)),))
-    end
-    return scalar(out)
+    codomain = inds(t; plev = 0)
+    return scalar(t * one(t, codomain, prime.(codomain)))
 end
 
 # One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through `project_aux`
@@ -518,43 +520,6 @@ end
 # `TensorAlgebra` generic, extended by ITensorBase for named tensors.
 using TensorAlgebra: matricize
 
-#
-# Diagonal manipulation. Legacy `map_diag(f, T)` applies `f` to the diagonal of a
-# (diagonal-like) tensor; used on factorization spectra (singular values /
-# eigenvalues). Reproduced via `MAK.diagview`, which aliases the diagonal storage on
-# every backend (a plain view for dense, blockwise for graded, per-sector for a
-# `TensorMap`), so no scalar indexing is needed.
-function _map_diagview!(f, dv)
-    copyto!(dv, map(f, dv))
-    return dv
-end
-# A backend's `diagview` may be a dict of per-block diagonal views (e.g. a plain
-# `TensorMap`, keyed by sector) rather than a single vector view.
-function _map_diagview!(f, dv::AbstractDict)
-    foreach(v -> _map_diagview!(f, v), values(dv))
-    return dv
-end
-function map_diag(f, T::AbstractITensor)
-    arr = copy(unnamed(T))
-    _map_diagview!(f, MAK.diagview(arr))
-    return nameddims(arr, ITensorBase.dimnames(T))
-end
-function map_diag!(f, T::AbstractITensor)
-    _map_diagview!(f, MAK.diagview(unnamed(T)))
-    return T
-end
-# Out-of-place-into-`dest` form `map_diag!(f, dest, src)`: write `f` of `src`'s diagonal
-# onto `dest`'s diagonal (TNQS calls it with `dest === src` for an in-place diagonal map).
-function map_diag!(f, dest::AbstractITensor, src::AbstractITensor)
-    d, s = unnamed(dest), unnamed(src)
-    if d === s
-        map_diag!(f, dest)
-    else
-        copyto!(MAK.diagview(d), map(f, MAK.diagview(s)))
-    end
-    return dest
-end
-
 #
 # Storage / element type accessors.
 #
diff --git a/src/imports.jl b/src/imports.jl
index 04ba865..adf8961 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -56,9 +56,8 @@ import .ITensorsITensorBaseCompat as ITensors
 using .ITensorsITensorBaseCompat:
     inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, cat_inds,
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
-    itensor, random_itensor, scalar, delta, id, onehot,
+    itensor, random_itensor, scalar, delta, similar_map, onehot,
     qr, svd, eigen, factorize, factorize_svd,
-    map_diag, map_diag!,
     scalartype, datatype, array, data,
     denseblocks, dense, hasqns,
     contract, inner, apply,
diff --git a/src/rdm.jl b/src/rdm.jl
index 8a40045..156a822 100644
--- a/src/rdm.jl
+++ b/src/rdm.jl
@@ -1,8 +1,6 @@
 function normalize_rdm(ρ::ITensor)
-    tr_ρ = copy(ρ)
-    for i in inds(ρ; plev = 0)
-        tr_ρ *= id(scalartype(ρ), (dag(i),), (dag(prime(i)),))
-    end
+    codomain = inds(ρ; plev = 0)
+    tr_ρ = ρ * one(ρ, codomain, prime.(codomain))
     return ρ / scalar(tr_ρ)
 end
 

From 379e83a0d801af0c66f63347609fcd0e39c36c32 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 12:05:37 -0400
Subject: [PATCH 08/14] Use the compatibility layer for op in the
 Heisenberg-picture example

The example loaded the legacy `ITensors` package for `op` and `ITensor`. Bind them to the compatibility layer through `TensorNetworkQuantumSimulator` instead, matching the thermal-state example, so no example pulls in legacy `ITensors`.
---
 examples/2dIsing_dynamics_Heisenbergpicture.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/2dIsing_dynamics_Heisenbergpicture.jl b/examples/2dIsing_dynamics_Heisenbergpicture.jl
index f9e6f51..3f1801d 100644
--- a/examples/2dIsing_dynamics_Heisenbergpicture.jl
+++ b/examples/2dIsing_dynamics_Heisenbergpicture.jl
@@ -1,7 +1,7 @@
 using TensorNetworkQuantumSimulator
 using Graphs: center
 using TensorNetworkQuantumSimulator: setindex_preserve!, noprime
-using ITensors: ITensors, ITensor
+using TensorNetworkQuantumSimulator: ITensors, ITensor
 
 function main()
     nx, ny = 4, 4

From 9d6a6c69ee0be0445b0633919646f929318b089b Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 12:30:26 -0400
Subject: [PATCH 09/14] Carry the dual arrows in the device-following identity
 contractions

The trace identity (`tr`, `normalize_rdm`) and the bilinear-form operator now build over the duals of the contracted legs, `dag.(codomain)` and `dag.(prime.(codomain))`, through `similar_map`. The device-following identity migration had dropped the `dag` that the previous construction carried. Dense backends ignore arrows, so that migration's tests passed, but a graded contraction requires arrow-opposite axes, so `reduced_density_matrix`, `tr`, and `inner` failed with an axis mismatch. The fix stays device-following and contracts on both dense and graded backends.
---
 src/Forms/bilinearform.jl                 |  6 +++++-
 src/ITensorsITensorBaseCompat/itensors.jl | 14 ++++++++------
 src/rdm.jl                                |  5 +++--
 3 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/src/Forms/bilinearform.jl b/src/Forms/bilinearform.jl
index aa199d5..faadbdb 100644
--- a/src/Forms/bilinearform.jl
+++ b/src/Forms/bilinearform.jl
@@ -18,7 +18,11 @@ function BilinearForm(ket::TensorNetworkState, bra::TensorNetworkState)
     sinds = siteinds(ket)
     verts = collect(vertices(ket))
     bra = TensorNetworkState(Dictionary(verts, [bra_tensor(bra, v) for v in verts]))
-    operator_tensors = [one(similar_map(ket[v], sinds[v], prime.(sinds[v])), sinds[v], prime.(sinds[v])) for v in verts]
+    operator_tensors = [
+        let codomain = dag.(sinds[v]), domain = dag.(prime.(sinds[v]))
+            one(similar_map(ket[v], codomain, domain), codomain, domain)
+        end for v in verts
+    ]
     operator = TensorNetworkState(Dictionary(verts, operator_tensors))
     return BilinearForm(ket, operator, bra)
 end
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 4d4d0a3..1ed4552 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -264,13 +264,15 @@ delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is))
 # pairs every unprimed index with its prime, the same construction `normalize_rdm` uses. This
 # is the index-paired definition — independent of storage order — rather than `sum` of the
 # underlying array's storage diagonal, which is only the trace for a rank-2 tensor whose legs
-# happen to be ordered to align. The domain is `prime.(codomain)` (not `inds(t; plev=1)`) so
-# each codomain index pairs with its own prime; `plev` filtering does not preserve that
-# pairing order. Accessed qualified (`ITensors.tr`) so it doesn't shadow `LinearAlgebra.tr`,
-# which TNQS still calls on plain matrices.
+# happen to be ordered to align. The identity carries the duals of `t`'s legs (`dag`) so each
+# leg contracts its partner on a graded backend, and its domain is `dag.(prime.(codomain))`
+# (not `inds(t; plev=1)`) so each codomain index pairs with its own prime; `plev` filtering
+# does not preserve that pairing order. Accessed qualified (`ITensors.tr`) so it doesn't shadow
+# `LinearAlgebra.tr`, which TNQS still calls on plain matrices.
 function tr(t::AbstractITensor)
-    codomain = inds(t; plev = 0)
-    return scalar(t * one(t, codomain, prime.(codomain)))
+    unprimed = inds(t; plev = 0)
+    codomain, domain = dag.(unprimed), dag.(prime.(unprimed))
+    return scalar(t * one(similar_map(t, codomain, domain), codomain, domain))
 end
 
 # One-hot vector along `i` at position `p` (legacy `onehot(i => p)`), through `project_aux`
diff --git a/src/rdm.jl b/src/rdm.jl
index 156a822..ff47950 100644
--- a/src/rdm.jl
+++ b/src/rdm.jl
@@ -1,6 +1,7 @@
 function normalize_rdm(ρ::ITensor)
-    codomain = inds(ρ; plev = 0)
-    tr_ρ = ρ * one(ρ, codomain, prime.(codomain))
+    unprimed = inds(ρ; plev = 0)
+    codomain, domain = dag.(unprimed), dag.(prime.(unprimed))
+    tr_ρ = ρ * one(similar_map(ρ, codomain, domain), codomain, domain)
     return ρ / scalar(tr_ρ)
 end
 

From dc5f2ccdbf3ea8b6b4a60f265882e85949da63ae Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 13:30:39 -0400
Subject: [PATCH 10/14] Match gate site indices by name so apply_gates works on
 graded backends

`vertices(::ITensor, ::TensorNetworkState)` resolved a gate's vertices with a full-`Index` `intersect`, which is arrow-sensitive on a graded axis: a gate's contraction leg is the dual of the site index, so the match came up empty and `apply_gates` hit a `BoundsError` in `simple_update`. The dense backend passed only because it ignores arrows. It now matches by `IndexName` through `nameisdisjoint`.

This also reworks the compat index-set algebra into a generic small-collection layer keyed by a `by` transform (`smallintersect` and friends, a linear scan with no hashing), name-keyed wrappers built on it with `by = name`, and the `commoninds`/etc. named-tensor family delegating to those.
---
 src/ITensorsITensorBaseCompat/itensors.jl | 61 +++++++++++++++++------
 src/TensorNetworks/tensornetworkstate.jl  |  5 +-
 2 files changed, 47 insertions(+), 19 deletions(-)

diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 1ed4552..356dc57 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -59,22 +59,51 @@ _compat_inds(t::AbstractITensor) = inds(t)
 _compat_inds(is) = is
 
 #
-# Index-set algebra. Legacy ITensors provides these as set operations over the
-# indices of its arguments; in the next-gen stack the same falls out of Base set ops
-# on `inds` (index identity/equality is by id).
+# Small-collection set operations, keyed by a transform `by` (like `sort`/`unique`'s `by`):
+# elements compare equal when `by(x) == by(y)`. Base's `intersect`/`union`/`setdiff`/… build
+# `Set`s (hashing); these instead scan `by(x) ∈ Iterators.map(by, b)`, a linear
+# O(length(a)*length(b)) pass over a lazy non-allocating view of the keys, the right tradeoff
+# when the collections are a handful of elements. The intersect/setdiff/union/symdiff forms
+# return elements of the first argument, always as `Vector`s.
 #
-# Matching is by index *name* (`IndexName`), not by full `Index` equality: two indices
-# name the same leg when their names match. On a graded axis a shared bond appears as an
-# index on one tensor and its dual (`conj`) on the other — same name, opposite arrow — so
-# `==` on the full `Index` misses it while the names still match. On the dense backend
-# name-matching and `==` coincide, so this is a strict generalization. The full `Index` is
-# returned, not the bare name: callers take `dim`/`tags` of the result and feed it to
-# `Index`/`qr`/`replaceinds`, which need the space. Results are always `Vector`s.
-_indvec(x) = collect(_compat_inds(x))
-commoninds(a, b) = (nb = name.(_indvec(b)); filter(i -> name(i) ∈ nb, _indvec(a)))
-uniqueinds(a, b) = (nb = name.(_indvec(b)); filter(i -> name(i) ∉ nb, _indvec(a)))
-unioninds(a, b) = vcat(_indvec(a), uniqueinds(b, a))
-hascommoninds(a, b) = !isempty(commoninds(a, b))
+smallintersect(a, b; by = identity) = (kb = Iterators.map(by, b); [x for x in a if by(x) ∈ kb])
+smallsetdiff(a, b; by = identity) = (kb = Iterators.map(by, b); [x for x in a if by(x) ∉ kb])
+smallunion(a, b; by = identity) = vcat(collect(a), smallsetdiff(b, a; by))
+smallsymdiff(a, b; by = identity) = vcat(smallsetdiff(a, b; by), smallsetdiff(b, a; by))
+smallisdisjoint(a, b; by = identity) = (kb = Iterators.map(by, b); !any(x -> by(x) ∈ kb, a))
+smallissubset(a, b; by = identity) = (kb = Iterators.map(by, b); all(x -> by(x) ∈ kb, a))
+smallissetequal(a, b; by = identity) = smallissubset(a, b; by) && smallissubset(b, a; by)
+
+#
+# Name-based index-set algebra: the small-collection ops keyed by `IndexName` (`by = name`)
+# rather than full `Index` equality. On a graded axis a shared bond appears as an index on one
+# tensor and its dual (`conj`) on the other, same name and opposite arrow, so `==` on the full
+# `Index` misses it while the names match; on the dense backend name-matching and `==` coincide,
+# so this is a strict generalization. The `Index`-returning ops give back the full `Index`es
+# (callers need the space to feed `Index`/`qr`/`replaceinds`). This is the shape we want upstream
+# in ITensorBase.
+#
+# Pairwise name equality (two single indices name the same leg).
+nameisequal(i, j) = name(i) == name(j)
+nameintersect(a, b) = smallintersect(a, b; by = name)
+namesetdiff(a, b) = smallsetdiff(a, b; by = name)
+nameunion(a, b) = smallunion(a, b; by = name)
+namesymdiff(a, b) = smallsymdiff(a, b; by = name)
+nameisdisjoint(a, b) = smallisdisjoint(a, b; by = name)
+nameissubset(a, b) = smallissubset(a, b; by = name)
+nameissetequal(a, b) = smallissetequal(a, b; by = name)
+
+#
+# Named-tensor index-set algebra (the `commoninds`/etc. surface). `_compat_inds` coerces a
+# tensor to its indices and passes an index collection through unchanged, so these currently
+# accept either; the `name*` ops take it from there (and always return `Vector`s, so a
+# returned index set can be fed straight back in, as `reduce(noncommoninds, tensors)` does).
+# These are meant to go back upstream to ITensorBase with stricter tensor-only definitions.
+#
+commoninds(a, b) = nameintersect(_compat_inds(a), _compat_inds(b))
+uniqueinds(a, b) = namesetdiff(_compat_inds(a), _compat_inds(b))
+unioninds(a, b) = nameunion(_compat_inds(a), _compat_inds(b))
+hascommoninds(a, b) = !nameisdisjoint(_compat_inds(a), _compat_inds(b))
 
 # Singular forms: the one common / one unique index. `noncommonind` is used in TNQS
 # as "the index of `a` not shared with `b`" (e.g. the non-contracted leg of an
@@ -82,7 +111,7 @@ hascommoninds(a, b) = !isempty(commoninds(a, b))
 commonind(a, b) = (cs = commoninds(a, b); isempty(cs) ? nothing : first(cs))
 noncommonind(a, b) = (us = uniqueinds(a, b); isempty(us) ? nothing : first(us))
 # Plural: indices not shared by both (symmetric difference).
-noncommoninds(a, b) = vcat(uniqueinds(a, b), uniqueinds(b, a))
+noncommoninds(a, b) = namesymdiff(_compat_inds(a), _compat_inds(b))
 
 # Index dimension (legacy `dim`). `dim(i)` is the length; `dim(is)` the product.
 dim(i::Index) = length(i)
diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl
index b4140ed..ad19029 100644
--- a/src/TensorNetworks/tensornetworkstate.jl
+++ b/src/TensorNetworks/tensornetworkstate.jl
@@ -1,4 +1,4 @@
-using .ITensorsITensorBaseCompat: random_itensor
+using .ITensorsITensorBaseCompat: nameisdisjoint, random_itensor
 
 """
     TensorNetworkState{V} <: AbstractTensorNetwork{V}
@@ -220,6 +220,5 @@ function tensornetworkstate(f::Function, args...)
 end
 
 function NamedGraphs.vertices(t::ITensor, tns::TensorNetworkState)
-    t_inds = inds(t)
-    return filter(v -> !isempty(intersect(t_inds, siteinds(tns, v))), collect(vertices(tns)))
+    return filter(v -> !nameisdisjoint(inds(t), siteinds(tns, v)), collect(vertices(tns)))
 end

From 8d1691bb1617e5df6ee7665f0160d3097297d3ab Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 14:13:32 -0400
Subject: [PATCH 11/14] Build order-2 identities via id/one and drop the
 order-2 delta method

`op("I")` and the form BP `default_message` built their order-2 identity maps through `delta`/`project`. They now use the identity-map primitives instead: `op("I")` uses the from-scratch `id(eltype, codomain, domain)` (only the site index is in hand, no prototype tensor to follow), and `default_message(form)` uses the device-following `one(similar_map(...), ...)` pattern over the ket bonds (codomain) and the duals of the bra bonds (domain), which generalizes to any number of ket-to-bra bond pairs. With no order-2 identity caller left on `delta`, the `delta(eltype, ::Tuple{Index,Index})` method is deleted. The remaining `delta` callers are dense copy tensors that fall to the general `diagonaltensor` path.
---
 src/Forms/abstractform.jl                 |  4 +++-
 src/ITensorsITensorBaseCompat/itensors.jl | 24 +++++++++++------------
 src/ITensorsITensorBaseCompat/ops.jl      |  5 +----
 3 files changed, 16 insertions(+), 17 deletions(-)

diff --git a/src/Forms/abstractform.jl b/src/Forms/abstractform.jl
index e1c9ebd..70b4c2d 100644
--- a/src/Forms/abstractform.jl
+++ b/src/Forms/abstractform.jl
@@ -22,7 +22,9 @@ function virtualinds(form::AbstractForm, edge::NamedEdge)
 end
 
 function default_message(form::AbstractForm, edge::AbstractEdge)
-    return adapt_like(form, denseblocks(delta(virtualinds(form, edge))))
+    cod = virtualinds(ket(form), edge)
+    dom = dag.(bra_virtualinds(form, edge))
+    return one(similar_map(ket(form)[src(edge)], cod, dom), cod, dom)
 end
 
 function bp_factors(form::AbstractForm, verts::Vector)
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 356dc57..0b96d7c 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -270,19 +270,19 @@ function similar_map(prototype::AbstractITensor, codomain, domain)
     return similar_map(prototype, scalartype(prototype), codomain, domain)
 end
 
+# From-scratch identity map: a dense identity embedded onto the `codomain`/`domain` index
+# partition via checked `project`, so the index axes select the backend (dense, graded,
+# `TensorMap`). Unlike `one(a, codomain, domain)` it needs no prototype tensor, so it is the
+# right primitive when only the indices and an element type are in hand (e.g. `op("I")`).
+function id(eltype::Type, codomain, domain)
+    m = Matrix{eltype}(LinearAlgebra.I, prod(length, codomain), prod(length, domain))
+    return TensorAlgebra.project(m, Tuple(codomain), Tuple(domain))
+end
+
+# Dense Kronecker copy (`delta`) tensor over the index axes. Dense-only: a super-diagonal
+# generally cannot be embedded while preserving a nontrivial symmetry, so graded/`TensorMap`
+# callers that want an order-2 identity build it via `id`/`one` at the callsite instead.
 delta(eltype::Type, is::Tuple) = diagonaltensor(ones(eltype, minimum(length, is)), is)
-# The order-2 delta is the identity map over the index pair, built via checked `project` so
-# the index axes select the backend (dense, graded, `TensorMap`). Legacy QN-`delta`
-# convention: the pair is arrow-opposite and the caller passes `(i, j)` with `j` un-dualized.
-# It coincides with the order-2 identity map, but stays its own method (its own calling
-# convention, and the count dispatch below routes graded default messages here). Higher-order deltas stay
-# on the dense `diagonaltensor` path above (a super-diagonal generally cannot be embedded
-# while preserving a nontrivial symmetry).
-function delta(eltype::Type, is::Tuple{Index, Index})
-    i, j = is
-    m = Matrix{eltype}(LinearAlgebra.I, length(i), length(j))
-    return TensorAlgebra.project(m, (i,), (dag(j),))
-end
 delta(eltype::Type, is::Index...) = delta(eltype, is)
 delta(eltype::Type, is::AbstractVector{<:Index}) = delta(eltype, Tuple(is))
 delta(is::Tuple) = delta(Float64, is)
diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/ITensorsITensorBaseCompat/ops.jl
index 14ee995..ed80fdb 100644
--- a/src/ITensorsITensorBaseCompat/ops.jl
+++ b/src/ITensorsITensorBaseCompat/ops.jl
@@ -94,10 +94,7 @@ end
 function op(name::AbstractString, sites::Index...; kwargs...)
     if name in ("I", "Id") && length(sites) == 1
         s = only(sites)
-        return project(
-            Matrix{ComplexF64}(LinearAlgebra.I, length(s), length(s)),
-            (ITensorBase.prime(s),), (s,)
-        )
+        return id(ComplexF64, (ITensorBase.prime(s),), (s,))
     end
     return op(OpName(name), SiteType("S=1/2"), sites...; kwargs...)
 end

From 01f712d8976fe32d1461601c7ccce05a592b2212 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 15:33:56 -0400
Subject: [PATCH 12/14] Use name-based namesetdiff for index set differences in
 qr, auxinds, and full_update

qr, auxinds, and full_update computed index complements with Base setdiff, which
builds a Set of Indexes. On a hash collision Julia falls back to
isequal(::Index, ::Index), which iterates the index space, and a TensorKit
GradedSpace is not iterable, so these threw on symmetric backends once a collision
occurred (small index collections happened to work, larger ones did not). They now
use the name-based namesetdiff, matching the rest of the compat name-keyed
index-set algebra, which never hashes or iterates an Index.
---
 src/Apply/full_update.jl                  | 5 +++--
 src/ITensorsITensorBaseCompat/itensors.jl | 2 +-
 src/TensorNetworks/tensornetworkstate.jl  | 4 ++--
 3 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl
index e2f626c..08edb76 100644
--- a/src/Apply/full_update.jl
+++ b/src/Apply/full_update.jl
@@ -1,4 +1,5 @@
 using KrylovKit: linsolve
+using .ITensorsITensorBaseCompat: namesetdiff
 
 """
     full_update(o::ITensor, ψ::TensorNetworkState, v⃗; envs, kwargs...)
@@ -115,8 +116,8 @@ function optimise_p_q(
 
     fstart = print_fidelity_loss ? fidelity(envs, p_cur, q_cur, p, q, o) : 0
 
-    qs_ind = setdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur)))))
-    ps_ind = setdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur)))))
+    qs_ind = namesetdiff(inds(q_cur), collect(Iterators.flatten(inds.(vcat(envs, p_cur)))))
+    ps_ind = namesetdiff(inds(p_cur), collect(Iterators.flatten(inds.(vcat(envs, q_cur)))))
 
     function b(p::ITensor, q::ITensor, o::ITensor, envs::Vector{ITensor}, r::ITensor)
         ts = vcat(ITensor[p, q, o, dag(prime(r))], envs)
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index 0b96d7c..a911100 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -327,7 +327,7 @@ const svd_trunc = MAK.svd_trunc
 # `linds` may be splatted, a single index, or one collection.
 function qr(a::AbstractITensor, linds...)
     left = cat_inds(linds...)
-    right = setdiff(collect(inds(a)), left)
+    right = namesetdiff(inds(a), left)
     return MAK.qr_compact(a, Tuple(left), Tuple(right))
 end
 
diff --git a/src/TensorNetworks/tensornetworkstate.jl b/src/TensorNetworks/tensornetworkstate.jl
index ad19029..01ee66b 100644
--- a/src/TensorNetworks/tensornetworkstate.jl
+++ b/src/TensorNetworks/tensornetworkstate.jl
@@ -1,4 +1,4 @@
-using .ITensorsITensorBaseCompat: nameisdisjoint, random_itensor
+using .ITensorsITensorBaseCompat: nameisdisjoint, namesetdiff, random_itensor
 
 """
     TensorNetworkState{V} <: AbstractTensorNetwork{V}
@@ -65,7 +65,7 @@ end
 bra_tensor(tns::TensorNetworkState, v) = bra_tensor(tns[v], auxinds(tns, v))
 
 # The dangling non-physical legs of a vertex tensor: dangling legs that are not site indices.
-auxinds(tns::TensorNetworkState, v) = Index[i for i in setdiff(uniqueinds(tns, v), siteinds(tns, v))]
+auxinds(tns::TensorNetworkState, v) = Index[i for i in namesetdiff(uniqueinds(tns, v), siteinds(tns, v))]
 
 # `auxinds_f` overrides the live aux-leg classification. The loop-correction weights
 # need this: they deliberately relabel a bond so it dangles in the modified network,

From a13dfde3f12f128e1119d9cb23f31c489e468164 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 18:22:12 -0400
Subject: [PATCH 13/14] Inline balanced SVD in the apply routines and simplify
 compat truncation

Replaces the multi-purpose compat `factorize_svd` with a direct balanced SVD inlined at the two apply call sites. `simple_update` and `full_update` now call `svd_trunc` and absorb the square root of the singular values into both factors, reading off the singular values and truncation error inline. That is all `factorize_svd` did for its only two callers, which both used the balanced (non-orthogonal) split.

Replaces `_trunc_spec` with a plain `itensor_trunc(; maxdim, cutoff)` that builds a MatrixAlgebraKit truncation strategy (`maxdim` caps the rank, `cutoff` maps to a relative 2-norm error) and returns `notrunc()` when neither applies. `generic_apply`'s no-cap default becomes `nothing` rather than `typemax(Int)`, so no sentinel handling leaks into the utility, and the orthogonalizing `factorize` wrapper uses `itensor_trunc` too.

Folds in a few smaller compat cleanups: `normalize_rdm` goes through `ITensors.tr`, the reshaping `itensor` constructor throws `DimensionMismatch` on a length mismatch, and the internal hermitian eigen helper is renamed and relies on `eigh_full`'s own hermitian tolerance instead of pre-projecting onto the hermitian part. Also trims over-long and outdated comments in the compat layer.

The truncation error is still computed from norms, which loses precision when little weight is discarded. Surfacing MatrixAlgebraKit's own truncation error through `ITensorBase.svd_trunc` is left for a follow-up, noted at both call sites.
---
 src/Apply/full_update.jl                      |  23 +-
 src/Apply/simple_update.jl                    |  26 +-
 .../ITensorsITensorBaseCompat.jl              |   2 +-
 src/ITensorsITensorBaseCompat/itensors.jl     | 260 +++++-------------
 src/ITensorsITensorBaseCompat/ops.jl          |  26 +-
 src/MessagePassing/boundarympscache.jl        |   2 +-
 src/imports.jl                                |   2 +-
 src/rdm.jl                                    |   5 +-
 8 files changed, 118 insertions(+), 228 deletions(-)

diff --git a/src/Apply/full_update.jl b/src/Apply/full_update.jl
index 08edb76..75b81fc 100644
--- a/src/Apply/full_update.jl
+++ b/src/Apply/full_update.jl
@@ -38,15 +38,20 @@ function full_update(
         apply_kwargs...,
     )
     if symmetrize
-        singular_values! = Ref{ITensor}()
-        Rᵥ₁, Rᵥ₂, spec = factorize_svd(
-            Rᵥ₁ * Rᵥ₂,
-            inds(Rᵥ₁);
-            ortho = "none",
-            singular_values!,
-            apply_kwargs...,
-        )
-        callback(; singular_values = singular_values![], truncation_error = spec.truncerr)
+        M = Rᵥ₁ * Rᵥ₂
+        codomain = inds(Rᵥ₁)
+        # Balanced SVD: split the singular values symmetrically (√S into each factor).
+        U, S, V = svd_trunc(M, codomain; trunc = itensor_trunc(; apply_kwargs...))
+        u = only(commoninds(U, S))
+        v = only(commoninds(S, V))
+        sqrtS = sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0)
+        Rᵥ₁, Rᵥ₂ = U * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V
+        # Best-effort truncation error from norms; suffers catastrophic cancellation when little is
+        # discarded. TODO: expose MatrixAlgebraKit's `ϵ` from `ITensorBase.svd_trunc` and use it here.
+        total = abs2(norm(M))
+        truncation_error = iszero(total) ? zero(real(scalartype(M))) :
+            max(zero(real(scalartype(M))), 1 - abs2(norm(S)) / total)
+        callback(; singular_values = S, truncation_error)
     end
     ψᵥ₁ = Qᵥ₁ * Rᵥ₁
     ψᵥ₂ = Qᵥ₂ * Rᵥ₂
diff --git a/src/Apply/simple_update.jl b/src/Apply/simple_update.jl
index 3fc5760..f4a3722 100644
--- a/src/Apply/simple_update.jl
+++ b/src/Apply/simple_update.jl
@@ -55,16 +55,22 @@ function simple_update(
         rᵥ₁ = commoninds(Qᵥ₁, Rᵥ₁)
         rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂)
         oR = ITensors.apply(o, Rᵥ₁ * Rᵥ₂)
-        singular_values! = Ref{ITensor}()
-        Rᵥ₁, Rᵥ₂, spec = factorize_svd(
-            oR,
-            unioninds(rᵥ₁, sᵥ₁);
-            ortho = "none",
-            singular_values!,
-            apply_kwargs...,
-        )
-        err = spec.truncerr
-        s_values = singular_values![]
+        # Balanced SVD: split the singular values symmetrically (√S into each factor) so neither
+        # side is isometric. The bond stays on `prime(u)` (keeping `u`'s name), so once this
+        # function `noprime`s its result the bond becomes `u`, which the returned `s_values` (over
+        # `(u, v)`) still shares for `apply_gate!`'s bond-message construction.
+        U, S, V = svd_trunc(oR, unioninds(rᵥ₁, sᵥ₁); trunc = itensor_trunc(; apply_kwargs...))
+        u = only(commoninds(U, S))
+        v = only(commoninds(S, V))
+        sqrtS = sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0)
+        Rᵥ₁, Rᵥ₂ = U * replaceind(sqrtS, v, prime(u)), replaceind(sqrtS, u, prime(u)) * V
+        s_values = S
+        # Best-effort truncation error from norms (SVD preserves the Frobenius norm); suffers
+        # catastrophic cancellation when little is discarded. TODO: expose MatrixAlgebraKit's `ϵ`
+        # from `ITensorBase.svd_trunc` and use it here instead.
+        total = abs2(norm(oR))
+        err = iszero(total) ? zero(real(scalartype(oR))) :
+            max(zero(real(scalartype(oR))), 1 - abs2(norm(S)) / total)
         Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)])
         Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)])
         updated_tensors = [Qᵥ₁ * Rᵥ₁, Qᵥ₂ * Rᵥ₂]
diff --git a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
index 21b475c..f38a706 100644
--- a/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
+++ b/src/ITensorsITensorBaseCompat/ITensorsITensorBaseCompat.jl
@@ -33,7 +33,7 @@ export
     # ITensor construction
     itensor, random_itensor, scalar, delta, onehot,
     # Factorizations
-    qr, svd, eigen, factorize, factorize_svd,
+    qr, svd, svd_trunc, eigen, factorize, itensor_trunc,
     # Storage / element-type accessors
     scalartype, datatype, array, data,
     # Dense / quantum-number no-ops
diff --git a/src/ITensorsITensorBaseCompat/itensors.jl b/src/ITensorsITensorBaseCompat/itensors.jl
index a911100..e2f7aa0 100644
--- a/src/ITensorsITensorBaseCompat/itensors.jl
+++ b/src/ITensorsITensorBaseCompat/itensors.jl
@@ -20,10 +20,8 @@ using LinearAlgebra: LinearAlgebra
 using TensorAlgebra.MatrixAlgebra: MatrixAlgebra
 using TensorAlgebra: TensorAlgebra
 
-# Legacy `inds(t; plev, tags)` accepted index-filtering keywords; the next-gen
-# `ITensorBase.inds` takes none. Own a compat `inds` that forwards to `ITensorBase.inds`
-# and applies the legacy filters — a compat-owned function (like `prime`/`dag`), not a
-# pirated method on `ITensorBase.inds`.
+# Legacy `inds(t; plev, tags)` took index-filtering keywords; `ITensorBase.inds` takes none. Compat
+# `inds` forwards to it and applies the legacy filters (a compat-owned function, not a pirated method).
 function inds(t::AbstractITensor; plev = nothing, tags = nothing)
     is = ITensorBase.inds(t)
     isnothing(plev) || (is = filter(i -> ITensorBase.plev(i) == plev, is))
@@ -31,20 +29,15 @@ function inds(t::AbstractITensor; plev = nothing, tags = nothing)
     return is
 end
 
-# Functions TNQS adds methods to for its own types (`scalartype`, `contract`,
-# `inner`, `datatype`). Rather than extend the upstream generics, TNQS owns these:
-# `scalartype` falls back to ITensorBase for tensors/arrays; `contract`/`inner` are
-# TNQS operations with their own base methods defined in the library. The per-file
-# method definitions extend these module-owned functions.
+# Module-owned functions TNQS extends for its own types: `scalartype` falls back to ITensorBase for
+# tensors/arrays; `contract`/`inner` are TNQS operations whose base methods live in the library.
 scalartype(x) = ITensorBase.scalartype(x)
 function inner end
 
-# Base contraction of a list of ITensors following a (possibly nested) pairwise
-# contraction `sequence` (legacy `ITensors.contract(tensors; sequence)`); leaves are
-# integer indices into `tensors`. TNQS adds the network-level `contract` methods.
-# The list is typed `AbstractVector` (not `AbstractVector{<:AbstractITensor}`) because
-# callers build it by concatenation that can widen to `Vector{Any}` (e.g. splicing in
-# an empty environment list).
+# Base contraction of a list of ITensors along a (possibly nested) pairwise `sequence` (legacy
+# `ITensors.contract(tensors; sequence)`); leaves are integer indices into `tensors`. Typed
+# `AbstractVector` (not `{<:AbstractITensor}`) because callers concatenate in ways that widen to
+# `Vector{Any}` (e.g. splicing in an empty environment list).
 function contract end
 function contract(tensors::AbstractVector; sequence = nothing)
     return isnothing(sequence) ? reduce(*, tensors) : _contract_sequence(tensors, sequence)
@@ -59,12 +52,11 @@ _compat_inds(t::AbstractITensor) = inds(t)
 _compat_inds(is) = is
 
 #
-# Small-collection set operations, keyed by a transform `by` (like `sort`/`unique`'s `by`):
-# elements compare equal when `by(x) == by(y)`. Base's `intersect`/`union`/`setdiff`/… build
-# `Set`s (hashing); these instead scan `by(x) ∈ Iterators.map(by, b)`, a linear
-# O(length(a)*length(b)) pass over a lazy non-allocating view of the keys, the right tradeoff
-# when the collections are a handful of elements. The intersect/setdiff/union/symdiff forms
-# return elements of the first argument, always as `Vector`s.
+# Small-collection set operations keyed by a transform `by` (elements compare equal when
+# `by(x) == by(y)`). These scan linearly rather than building `Set`s: Base's `Set`-based ops hash,
+# and hashing an `Index` over a TensorKit `GradedSpace` can fall back to iterating the space, which
+# throws. The intersect/setdiff/union/symdiff forms return elements of the first argument as
+# `Vector`s.
 #
 smallintersect(a, b; by = identity) = (kb = Iterators.map(by, b); [x for x in a if by(x) ∈ kb])
 smallsetdiff(a, b; by = identity) = (kb = Iterators.map(by, b); [x for x in a if by(x) ∉ kb])
@@ -75,15 +67,12 @@ smallissubset(a, b; by = identity) = (kb = Iterators.map(by, b); all(x -> by(x)
 smallissetequal(a, b; by = identity) = smallissubset(a, b; by) && smallissubset(b, a; by)
 
 #
-# Name-based index-set algebra: the small-collection ops keyed by `IndexName` (`by = name`)
-# rather than full `Index` equality. On a graded axis a shared bond appears as an index on one
-# tensor and its dual (`conj`) on the other, same name and opposite arrow, so `==` on the full
-# `Index` misses it while the names match; on the dense backend name-matching and `==` coincide,
-# so this is a strict generalization. The `Index`-returning ops give back the full `Index`es
-# (callers need the space to feed `Index`/`qr`/`replaceinds`). This is the shape we want upstream
-# in ITensorBase.
+# Name-based index-set algebra: the small-collection ops keyed by `IndexName` (`by = name`) rather
+# than full `Index` equality. On a graded axis a shared bond appears as an index on one tensor and
+# its dual (`conj`) on the other — same name, opposite arrow — so full-`Index` `==` misses it while
+# the names match; on the dense backend the two coincide. The `Index`-returning ops give back full
+# `Index`es (callers need the space to feed `Index`/`qr`/`replaceinds`).
 #
-# Pairwise name equality (two single indices name the same leg).
 nameisequal(i, j) = name(i) == name(j)
 nameintersect(a, b) = smallintersect(a, b; by = name)
 namesetdiff(a, b) = smallsetdiff(a, b; by = name)
@@ -94,11 +83,10 @@ nameissubset(a, b) = smallissubset(a, b; by = name)
 nameissetequal(a, b) = smallissetequal(a, b; by = name)
 
 #
-# Named-tensor index-set algebra (the `commoninds`/etc. surface). `_compat_inds` coerces a
-# tensor to its indices and passes an index collection through unchanged, so these currently
-# accept either; the `name*` ops take it from there (and always return `Vector`s, so a
-# returned index set can be fed straight back in, as `reduce(noncommoninds, tensors)` does).
-# These are meant to go back upstream to ITensorBase with stricter tensor-only definitions.
+# Named-tensor index-set algebra (the `commoninds`/etc. surface). `_compat_inds` coerces a tensor to
+# its indices and passes an index collection through unchanged, so these accept either; the `name*`
+# ops take it from there (always returning `Vector`s, so a result can be fed back in, as
+# `reduce(noncommoninds, tensors)` does).
 #
 commoninds(a, b) = nameintersect(_compat_inds(a), _compat_inds(b))
 uniqueinds(a, b) = namesetdiff(_compat_inds(a), _compat_inds(b))
@@ -157,20 +145,15 @@ replaceind(t, from::Index, to::Index) = replaceinds(t, from => to)
 _as_index_vec(x::Index) = [x]
 _as_index_vec(xs) = collect(xs)
 cat_inds(xs...) = reduce(vcat, map(_as_index_vec, xs))
-# Legacy `replaceinds` accepted collection arguments (`replaceinds(t, [i, k], [j, l])`
-# and `replaceinds(t, [i, k] => [j, l])`); ITensorBase provides only the pair-splat form
-# (`replaceinds(t, i => j, k => l)`). Own a compat `replaceinds` that handles the
-# collection forms and forwards the pair-splat form.
+# Legacy `replaceinds` took collection arguments (`replaceinds(t, [i,k], [j,l])`, `... => ...`);
+# ITensorBase provides only the pair-splat form, so this compat handles the collection forms.
 #
-# The base case relabels *names* via `replacedimnames`, not `ITensorBase.replaceinds`.
-# The latter replaces the index *space* (routing through `axes`/`getindex`), which
-# scalar-indexes a graded axis and errors; a pure name relabel is space-agnostic and works
-# on every backend. Keys are stripped to their `IndexName` because `replacedimnames`
-# matches on `dimnames` and silently no-ops on a raw `Index` key.
+# The base case relabels *names* via `replacedimnames`, not `ITensorBase.replaceinds` (which replaces
+# the index *space*, scalar-indexing a graded axis and erroring). Keys are stripped to `IndexName`
+# because `replacedimnames` matches on `dimnames` and silently no-ops on a raw `Index` key.
 #
-# NB: an `Index` is itself an `AbstractVector` (it's a `NamedUnitRange`), so the
-# collection types are constrained to `AbstractVector{<:Index}` to avoid capturing a
-# bare `Index` (an `AbstractVector{Int}`) and iterating over its integer range.
+# NB: an `Index` is itself an `AbstractVector` (`NamedUnitRange`), so the collection types are
+# constrained to `AbstractVector{<:Index}` to avoid capturing a bare `Index` and iterating its range.
 const _IndexColl = Union{Tuple{Vararg{Index}}, AbstractVector{<:Index}}
 function replaceinds(t, pairs::Pair...)
     return ITensorBase.replacedimnames(
@@ -196,18 +179,18 @@ end
 # dimensions (e.g. a `d^2 × d^2` two-site gate matrix over four site legs).
 function itensor(array, is...)
     length(array) == prod(length, is) ||
-        error(
-        "array with $(length(array)) elements cannot fill indices of dimensions $(length.(is))"
+        throw(
+        DimensionMismatch(
+            "array with $(length(array)) elements cannot fill indices of dimensions $(length.(is))"
+        )
     )
     return reshape(array, map(length, is))[is...]
 end
 itensor(array, is::Union{Tuple, AbstractVector}) = itensor(array, is...)
 
-# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): adds a rank-0
-# `ITensor(x::Number)` constructor, which ITensorBase deliberately omits and does not plan to
-# support. Legacy ITensors uses it as a multiplicative identity to seed a product accumulator
-# (`out = ITensor(1); out *= t; ...`). Kept here for now; the accumulator call sites get rewritten
-# to a different pattern later, retiring this method rather than upstreaming it.
+# TYPE PIRACY (temporary, compat-owned — NOT upstream): rank-0 `ITensor(x::Number)`, which
+# ITensorBase deliberately omits. Legacy ITensors uses it to seed a product accumulator
+# (`out = ITensor(1); out *= t`). Retired once those call sites are rewritten.
 ITensorBase.ITensor(x::Number) = nameddims(fill(x), ())
 
 # Random ITensor over the given indices (legacy `random_itensor`).
@@ -219,12 +202,9 @@ random_itensor(is::Union{Tuple, AbstractVector}) = random_itensor(Float64, is)
 # Rank-0 scalar extraction (legacy `scalar`).
 scalar(t::AbstractITensor) = t[]
 
-# Dense Kronecker delta tensor. ITensorBase deliberately omits the `delta` tensor
-# type that legacy ITensors had; this dense version is vendored from
-# ITensorNetworksNext's `ITensorNetworkGenerators/delta_network.jl` (a dense delta
-# defined on a graph generator there), copied because TNQS doesn't depend on
-# ITensorNetworksNext for this migration. A graded/sector-aware `delta` is a stack
-# gap (tracked separately).
+# Dense Kronecker delta tensor (legacy `delta`), vendored from ITensorNetworksNext's
+# `ITensorNetworkGenerators/delta_network.jl` (TNQS doesn't depend on it for this migration). A
+# graded/sector-aware `delta` is a stack gap (tracked separately).
 diaglength(a::AbstractArray) = minimum(size(a))
 function diagstride(a::AbstractArray)
     s = 1
@@ -254,12 +234,10 @@ function diagonaltensor(
     return nameddims(diagonaltensor(diag, unnamed.(is)), name.(is))
 end
 
-# Allocate an identity-map-shaped tensor with undefined data over the `codomain`/`domain`
-# indices, following `prototype`'s backend and element type with the domain axes dualized in
-# storage (the `similar_map` map convention). `prototype` only donates the backend, so unlike
-# `Base.one(a, codomain, domain)` (which needs `a` to already carry the map's legs) a caller
-# whose prototype lacks a leg still builds a device-following identity via
-# `one(similar_map(prototype, codomain, domain), codomain, domain)`.
+# Allocate an identity-map-shaped tensor over the `codomain`/`domain` indices, following
+# `prototype`'s backend/eltype (domain axes dualized in storage). Unlike `Base.one(a, codomain,
+# domain)` — which needs `a` to already carry the map's legs — this lets a caller whose prototype
+# lacks a leg still build a device-following identity via `one(similar_map(...), codomain, domain)`.
 function similar_map(prototype::AbstractITensor, eltype::Type, codomain, domain)
     raw = TensorAlgebra.similar_map(
         unnamed(prototype), eltype, unnamed.(codomain), unnamed.(domain)
@@ -289,15 +267,10 @@ delta(is::Tuple) = delta(Float64, is)
 delta(is::Index...) = delta(Float64, is)
 delta(is::AbstractVector{<:Index}) = delta(Float64, Tuple(is))
 
-# Trace of an ITensor over its prime pairs (legacy `tr`): contract with the identity map that
-# pairs every unprimed index with its prime, the same construction `normalize_rdm` uses. This
-# is the index-paired definition — independent of storage order — rather than `sum` of the
-# underlying array's storage diagonal, which is only the trace for a rank-2 tensor whose legs
-# happen to be ordered to align. The identity carries the duals of `t`'s legs (`dag`) so each
-# leg contracts its partner on a graded backend, and its domain is `dag.(prime.(codomain))`
-# (not `inds(t; plev=1)`) so each codomain index pairs with its own prime; `plev` filtering
-# does not preserve that pairing order. Accessed qualified (`ITensors.tr`) so it doesn't shadow
-# `LinearAlgebra.tr`, which TNQS still calls on plain matrices.
+# Trace over prime pairs (legacy `tr`): contract with the identity map pairing each unprimed index
+# with its prime. The domain is built as `dag.(prime.(codomain))`, not `inds(t; plev=1)` — `plev`
+# filtering does not preserve the pairing order between the plev-0 and plev-1 groups. Accessed
+# qualified (`ITensors.tr`) so it doesn't shadow `LinearAlgebra.tr`, which TNQS calls on matrices.
 function tr(t::AbstractITensor)
     unprimed = inds(t; plev = 0)
     codomain, domain = dag.(unprimed), dag.(prime.(unprimed))
@@ -358,7 +331,7 @@ end
 # truncation kwargs (`cutoff`/`maxdim`/`mindim`) are not yet translated to MAK's
 # `trunc=(; ...)` spec — see api_migration_map.md; the current callsites pass
 # `cutoff = nothing` (full decomposition).
-function _hermitian_eigh(
+function _eigh(
         m::AbstractITensor,
         codomain,
         domain;
@@ -371,12 +344,6 @@ function _hermitian_eigh(
         "the compat `eigen` does not yet translate the `cutoff` truncation kwarg to MatrixAlgebraKit's `trunc` spec"
     )
     cod, dom = _astuple(codomain), _astuple(domain)
-    # `MAK.eigh_full` rejects a matrix that is hermitian only up to numerical noise, but
-    # the caller asserts `ishermitian = true`. Project onto the hermitian part first,
-    # matching legacy `eigen`'s treat-as-hermitian behavior. `MAK.project_hermitian`
-    # works at the matricized level, so it stays generic over graded and `TensorMap`
-    # backends (a named-level `m + swap(conj(m))` is not).
-    m = MAK.project_hermitian(m, cod, dom)
     D, U = MAK.eigh_full(m, cod, dom)
     u = only(commoninds(D, U))         # eigenvalue index shared with U
     t = only(uniqueinds(D, U))         # D's other (independent) index
@@ -388,7 +355,7 @@ end
 # `m = Vt * D * dag(V)` (with `Vt` the relabeling of `V` from `Rinds` to `Linds`) —
 # no conjugation of `U`.
 function eigen(m::AbstractITensor, codomain, domain; kwargs...)
-    return _hermitian_eigh(m, codomain, domain; kwargs...)
+    return _eigh(m, codomain, domain; kwargs...)
 end
 
 # No-partition form `eigen(m)` matches legacy ITensors' `eigen(A)`, which auto-partitions
@@ -402,24 +369,20 @@ function eigen(m::AbstractITensor; kwargs...)
     length(p0) == length(p1) || error(
         "`eigen` without an index partition expects each plev-0 index to be paired with its prime"
     )
-    D, U = _hermitian_eigh(m, Tuple(p1), Tuple(p0); kwargs...)
+    D, U = _eigh(m, Tuple(p1), Tuple(p0); kwargs...)
     return D, conj(U)
 end
 
-# Translate the legacy `cutoff`/`maxdim` truncation kwargs to a MatrixAlgebraKit
-# `trunc` strategy. ITensors `cutoff` discards the smallest singular values whose
-# summed squares are a fraction ≤ cutoff of the total, i.e. a relative 2-norm
-# truncation error of `sqrt(cutoff)` on the singular-value vector — MAK's
-# `truncerror(; rtol = sqrt(cutoff), p = 2)`. `maxdim` caps the kept rank
-# (`truncrank`). When both apply, the intersection keeps the more aggressive of the
-# two, matching ITensors. Returns `nothing` when neither truncates.
-function _trunc_spec(cutoff, maxdim)
-    specs = []
-    isnothing(maxdim) || maxdim ≥ typemax(Int) || push!(specs, MAK.truncrank(maxdim))
-    isnothing(cutoff) || iszero(cutoff) ||
-        push!(specs, MAK.truncerror(; rtol = sqrt(cutoff), p = 2))
-    isempty(specs) && return nothing
-    return reduce(&, specs)
+# ITensors-style truncation as a MatrixAlgebraKit strategy. `maxdim` caps the kept rank
+# (`truncrank`); `cutoff` discards the smallest singular values whose summed squares are a fraction
+# ≤ cutoff of the total, i.e. a relative 2-norm error of `sqrt(cutoff)` on the singular-value vector
+# (`truncerror(; rtol = sqrt(cutoff), p = 2)`). Either kwarg may be `nothing` (or trivial) to drop
+# that constraint; with neither, the result is `notrunc()`.
+function itensor_trunc(; maxdim = nothing, cutoff = nothing)
+    trunc = MAK.notrunc()
+    isnothing(maxdim) || (trunc &= MAK.truncrank(maxdim))
+    isnothing(cutoff) || iszero(cutoff) || (trunc &= MAK.truncerror(; rtol = sqrt(cutoff), p = 2))
+    return trunc
 end
 
 # Split `a`'s indices into (left, right) by the requested `linds`, matching by name:
@@ -450,8 +413,8 @@ function factorize(
         tags = nothing
     )
     left, right = _bipartition_inds(a, linds)
-    trunc = _trunc_spec(cutoff, maxdim)
-    if isnothing(trunc)
+    notruncation = (isnothing(cutoff) || iszero(cutoff)) && isnothing(maxdim)
+    if notruncation
         if ortho == "left"
             L, R = MAK.qr_compact(a, Tuple(left), Tuple(right))
         elseif ortho == "right"
@@ -462,7 +425,7 @@ function factorize(
             )
         end
     else
-        U, S, Vt = MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc)
+        U, S, Vt = MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc = itensor_trunc(; cutoff, maxdim))
         if ortho == "left"
             L, R = U, S * Vt
         elseif ortho == "right"
@@ -481,67 +444,7 @@ function factorize(
     return L, R
 end
 
-# Absorb an SVD `(U, S, Vt)` into two factors per `ortho`: `"left"` makes `U` isometric,
-# `"right"` makes `Vt` isometric, `"none"` splits `S = √S · √S` so the weight is shared.
-# For `"none"` the bond lands on `prime(u)` (where `u`/`v` are the SVD's left/right
-# indices) and the separately-returned `singular_values!` stays over the unprimed
-# `(u, v)` — matching legacy ITensors so a caller that `noprime`s the factors (e.g.
-# `simple_update`) ends up with the bond on `u`, which `S` still shares.
-function _absorb_svd(U, S, Vt, ortho)
-    if ortho == "left"
-        return U, S * Vt
-    elseif ortho == "right"
-        return U * S, Vt
-    elseif ortho == "none"
-        u = only(commoninds(U, S))
-        v = only(commoninds(S, Vt))
-        up = prime(u)
-        # `S` is diagonal, so this hits the diagonal fast path (no eigendecomposition)
-        # on every backend. No clamping: zero singular values stay zero under `^(1//2)`.
-        sqrtS = MatrixAlgebra.sqrth_safe(S, (u,), (v,); atol = 0, rtol = 0)
-        return U * replaceind(sqrtS, v, up), Vt * replaceind(sqrtS, u, up)
-    end
-    return error(
-        "compat `factorize_svd` supports ortho = \"left\" / \"right\" / \"none\" (got $(repr(ortho)))"
-    )
-end
 
-# Legacy `factorize_svd(a, linds...; ortho, singular_values!, cutoff, maxdim, tags)`:
-# an always-SVD factorization returning `(L, R, spec)`. `spec.truncerr` is the fraction
-# of squared spectral weight discarded (SVD preserves the Frobenius norm, so the total
-# weight is `norm(a)^2`). If `singular_values!` is a `Ref`, it is filled with `S`.
-function factorize_svd(
-        a::AbstractITensor,
-        linds...;
-        ortho = "left",
-        singular_values! = nothing,
-        cutoff = nothing,
-        maxdim = nothing,
-        tags = nothing
-    )
-    left, right = _bipartition_inds(a, linds)
-    trunc = _trunc_spec(cutoff, maxdim)
-    U, S, Vt = if isnothing(trunc)
-        MAK.svd_compact(a, Tuple(left), Tuple(right))
-    else
-        MAK.svd_trunc(a, Tuple(left), Tuple(right); trunc)
-    end
-    total = abs2(LinearAlgebra.norm(a))
-    kept = abs2(LinearAlgebra.norm(S))
-    truncerr = if iszero(total)
-        zero(real(scalartype(a)))
-    else
-        max(zero(kept / total), 1 - kept / total)
-    end
-    isnothing(singular_values!) || (singular_values![] = S)
-    L, R = _absorb_svd(U, S, Vt, ortho)
-    if !isnothing(tags)
-        b = only(commoninds(L, R))
-        bnew = settags(b, tags)
-        L, R = replaceind(L, b, bnew), replaceind(R, b, bnew)
-    end
-    return L, R, (; truncerr)
-end
 
 #
 # Index fusion. The legacy `combiner` is retired (no compat shim); call sites fuse
@@ -563,17 +466,11 @@ datatype(T::AbstractITensor) = typeof(unnamed(T))
 array(T::AbstractITensor) = convert(Array, unnamed(T))
 data(T::AbstractITensor) = unnamed(T)
 
-# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): extends
-# `Adapt.adapt_structure` for `AbstractITensor` with an eltype target. Using
-# `Adapt.adapt_structure` for eltype *conversion* is an abuse of Adapt.jl (Adapt is for
-# storage/device adaptation, not changing the scalar type), so this does not belong upstream.
-# Kept here for now; the eltype-conversion call sites get rewritten with a different pattern
-# later, retiring this shim rather than upstreaming it.
-#
-# Legacy `adapt(eltype)(t)` converts an ITensor's scalar (element) type. ITensorBase's
-# Adapt integration adapts the storage array/device type but leaves the element type
-# alone, so reproduce the eltype conversion for a `Number` target (TNQS uses
-# `adapt(eltype)(state(...))` to build typed product states).
+# TYPE PIRACY (temporary, compat-owned — NOT upstream): extends `Adapt.adapt_structure` for
+# `AbstractITensor` with a `Number` eltype target, reproducing legacy `adapt(eltype)(t)` scalar-type
+# conversion (ITensorBase's Adapt integration adapts storage/device but leaves the eltype alone).
+# Using Adapt for eltype conversion is an abuse, so it does not belong upstream; retired once the
+# call sites move to a different pattern.
 function Adapt.adapt_structure(::Type{elt}, T::AbstractITensor) where {elt <: Number}
     # Short-circuit the no-op case: a non-`AbstractArray` backend (e.g. a `TensorMap`)
     # has no `convert(AbstractArray{elt}, ...)` method, but needs none when the element
@@ -601,10 +498,8 @@ hasqns(::Any) = false
 # The operator / named-state system (`op` / `state`) is vendored separately in
 # `ops.jl`, included right after this file by the module file.
 
-# Direct sum (legacy `directsum`): block-diagonal placement of several tensors along
-# specified axes, with the non-summed (shared) axes preserved. Vendored densely
-# because the next-gen stack has no `directsum` yet — a real missing ITensorBase
-# feature (tracked); drop this once it lands upstream.
+# Direct sum (legacy `directsum`): block-diagonal placement of tensors along the summed axes, shared
+# axes preserved. Vendored densely — the next-gen stack has no `directsum` yet (tracked upstream).
 #   directsum(out_inds, (t1 => summed_inds1), (t2 => summed_inds2), ...)
 function directsum(out_inds, pairs::Pair...)
     out_inds = Tuple(out_inds)
@@ -677,16 +572,11 @@ function settags(i::Index, d::AbstractDict)
     end
     return i
 end
-# TYPE PIRACY (temporary, compat-owned — NOT an upstream candidate): the two methods below add
-# `ITensorBase.Index` constructors taking a tag string / tag dict, a legacy positional-tag form
-# ITensorBase does not plan to support (the next-gen spelling passes tags via the `tags` keyword
-# argument). Kept here for now; the call sites get modernized to the `tags` kwarg later, retiring
-# these methods rather than upstreaming them.
-#
-# Legacy positional tagged-index constructor `Index(dim, "tag")`. Only the dimension
-# form is pirated: a range/space first argument collides with ITensorBase's own
-# two-argument `Index` constructors, so tagging an index minted over a backend axis
-# spells the two steps (`settags(Index(r), "tag")`).
+# TYPE PIRACY (temporary, compat-owned — NOT upstream): the two methods below add `Index`
+# constructors taking a tag string / tag dict (legacy positional-tag form; next-gen passes tags via
+# the `tags` kwarg). Retired once the call sites move to the kwarg form. Only the `dim` first-argument
+# form is pirated: a range/space first argument collides with ITensorBase's own two-argument `Index`
+# constructors, so tagging an axis-minted index spells the two steps (`settags(Index(r), "tag")`).
 ITensorBase.Index(dim::Integer, tagstr::AbstractString) = settags(Index(dim), tagstr)
 # Build a fresh index carrying a tag dictionary (legacy `Index(dim, tags(i))`, where
 # the next-gen `tags` returns a `Dict{String, String}`).
diff --git a/src/ITensorsITensorBaseCompat/ops.jl b/src/ITensorsITensorBaseCompat/ops.jl
index ed80fdb..0a17db5 100644
--- a/src/ITensorsITensorBaseCompat/ops.jl
+++ b/src/ITensorsITensorBaseCompat/ops.jl
@@ -1,14 +1,9 @@
-# Vendored minimal operator / named-state system (legacy `ITensors.op` / `ITensors.state`).
+# Vendored minimal operator / named-state system (legacy `ITensors.op` / `ITensors.state`), scoped to
+# the `S=1/2` single- and two-qubit gate set and the computational/Pauli-basis product states.
 #
-# Round-1 scope is the single-qubit (`S=1/2`) gate set and the computational/Pauli-basis
-# product states — enough for `siteinds` / `tensornetworkstate` state construction and
-# single-site `expect`. The full two-qubit gate set, the parametric and in-house gates,
-# and the Heisenberg/PTM (PauliPropagation) path belong to the gate-application round
-# (`Apply/`), which is still excluded; they are not vendored here.
-#
-# Legacy ITensors exposes these through the `OpName` / `SiteType` dispatch system. TNQS
-# only ever calls `op(name, sites...; kwargs...)` and `state(name, site)`, so the vendor
-# is a plain name-keyed lookup rather than a reimplementation of that dispatch machinery.
+# Legacy ITensors exposes these through the `OpName` / `SiteType` dispatch system. TNQS only ever
+# calls `op(name, sites...; kwargs...)` and `state(name, site)`, so this vendors a name-keyed lookup
+# over that machinery rather than reimplementing it.
 
 using TensorAlgebra: TensorAlgebra, project, tryproject
 
@@ -154,13 +149,10 @@ function op(::OpName"SWAP", ::SiteType"S=1/2")
     return Float64[1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
 end
 
-# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor`, inferring the
-# prime-pair codomain/domain and forwarding to ITensorBase's matricization
-# `exp(a, dimnames_codomain, dimnames_domain)` (graded-capable). Legacy ITensors provided
-# `exp(::ITensor)`; gates defined as `exp` of a Hamiltonian operator rely on it (e.g. a user
-# `op(::OpName"MyZRot", ...) = exp(-im θ/2 * op("Z", s))`). To de-pirate: make this compat-owned
-# (an `exp` in this module, not a `Base.exp` method). Not an upstream candidate (the upstream
-# matricization `exp` is the target, not a `Base.exp(::ITensor)`).
+# TYPE PIRACY (temporary): extends `Base.exp` for an operator `ITensor`, inferring the prime-pair
+# codomain/domain and forwarding to ITensorBase's matricization `exp(a, codomain, domain)`
+# (graded-capable). Gates defined as `exp` of a Hamiltonian operator rely on it (e.g. a user
+# `op(::OpName"MyZRot", ...) = exp(-im θ/2 * op("Z", s))`). De-pirate by making it compat-owned.
 function Base.exp(t::AbstractITensor)
     p0 = filter(i -> ITensorBase.plev(i) == 0, collect(inds(t)))
     isempty(p0) && error("exp(::ITensor) expects indices paired as (i, prime(i))")
diff --git a/src/MessagePassing/boundarympscache.jl b/src/MessagePassing/boundarympscache.jl
index 17ea521..fcef5d3 100644
--- a/src/MessagePassing/boundarympscache.jl
+++ b/src/MessagePassing/boundarympscache.jl
@@ -408,7 +408,7 @@ function generic_apply(
         mps::Dictionary{Int, <:ITensor},
         right_inds::Vector{<:Vector{<:Index}};
         cutoff = 0.0,
-        maxdim = typemax(Int),
+        maxdim = nothing,
         normalize = true,
     )
     b = length(mpo)
diff --git a/src/imports.jl b/src/imports.jl
index adf8961..c9c85c6 100644
--- a/src/imports.jl
+++ b/src/imports.jl
@@ -57,7 +57,7 @@ using .ITensorsITensorBaseCompat:
     inds, commoninds, commonind, uniqueinds, noncommonind, noncommoninds, unioninds, hascommoninds, cat_inds,
     sim, dag, prime, noprime, replaceind, replaceinds, dim, swapind,
     itensor, random_itensor, scalar, delta, similar_map, onehot,
-    qr, svd, eigen, factorize, factorize_svd,
+    qr, svd, svd_trunc, eigen, factorize, itensor_trunc,
     scalartype, datatype, array, data,
     denseblocks, dense, hasqns,
     contract, inner, apply,
diff --git a/src/rdm.jl b/src/rdm.jl
index ff47950..43b8ac1 100644
--- a/src/rdm.jl
+++ b/src/rdm.jl
@@ -1,8 +1,5 @@
 function normalize_rdm(ρ::ITensor)
-    unprimed = inds(ρ; plev = 0)
-    codomain, domain = dag.(unprimed), dag.(prime.(unprimed))
-    tr_ρ = ρ * one(similar_map(ρ, codomain, domain), codomain, domain)
-    return ρ / scalar(tr_ρ)
+    return ρ / ITensors.tr(ρ)
 end
 
 """

From 9d581b32aab5c76611aad7eac3c3641ce9ecf1a9 Mon Sep 17 00:00:00 2001
From: Matthew Fishman <mfishman@flatironinstitute.org>
Date: Tue, 7 Jul 2026 18:38:30 -0400
Subject: [PATCH 14/14] Remove [sources] pins now that the upstreams have
 registered

ITensorBase 0.10.9 and TensorAlgebra 0.17.0 are registered, so resolve
against the released versions instead of the branch pins. The [compat]
floors already require these versions.
---
 Project.toml | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 9f01eb3..5349d2c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -52,7 +52,3 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "Random"]
-
-[sources]
-ITensorBase = {url = "https://github.com/ITensor/ITensorBase.jl", rev = "mf/tensoralgebra-0.17"}
-TensorAlgebra = {url = "https://github.com/ITensor/TensorAlgebra.jl", rev = "mf/project-trailing-axes"}