From f716e1a8c60630e21d2ed927c4944f3ce3ef597c Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 7 Oct 2026 08:48:55 -0400 Subject: [PATCH 1/5] Use the asymmetric Gram root for the BP simple-update gauge Gauges each bond by the asymmetric root of its message and its inverse, rather than the symmetric root, which keeps the truncation optimal for a non-Hermitian environment. The inverse factor is twisted so that it inverts the forward factor under contraction rather than under the bare matrix product. Co-Authored-By: Claude Opus 5 (1M context) --- Project.toml | 10 +++++-- src/apply/apply_operators.jl | 54 ++++++++++++++++++++++++++++-------- test/Project.toml | 4 +-- 3 files changed, 52 insertions(+), 16 deletions(-) diff --git a/Project.toml b/Project.toml index 5b863052..deff89b0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorNetworksNext" uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c" -version = "0.11.0" +version = "0.11.1" authors = ["ITensor developers and contributors"] [workspace] @@ -27,6 +27,10 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" [weakdeps] OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715" +[sources.ITensorBase] +rev = "mf/named-twist" +url = "https://github.com/ITensor/ITensorBase.jl" + [extensions] ITensorNetworksNextOMEinsumContractionOrdersExt = "OMEinsumContractionOrders" @@ -38,7 +42,7 @@ Combinatorics = "1" DataGraphs = "0.6" Dictionaries = "0.4.5" Graphs = "1.13.1" -ITensorBase = "0.16" +ITensorBase = "0.16.1" LinearAlgebra = "1.10" MacroTools = "0.5.16" MatrixAlgebraKit = "0.6" @@ -47,5 +51,5 @@ OMEinsumContractionOrders = "1" Random = "1.10" SimpleTraits = "0.9.5" SplitApplyCombine = "1.2.3" -TensorAlgebra = "0.23" +TensorAlgebra = "0.23.1" julia = "1.10" diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl index 0682df57..5eaf68b9 100644 --- a/src/apply/apply_operators.jl +++ b/src/apply/apply_operators.jl @@ -2,11 +2,45 @@ using .AlgorithmsInterfaceExtensions: AlgorithmsInterfaceExtensions as AIE using AlgorithmsInterface: AlgorithmsInterface as AI using Base: @kwdef using Graphs: dst, src, vertices -using ITensorBase: AbstractITensor, apply, names, operator, rename +using ITensorBase: AbstractITensor, AbstractNamedTensor, ITensor, Index, NamedTensor, + NamedTensorOperator, apply, inputinds, inputnames, name, names, operator, outputnames, + rename, sim, state, to_inds, uniquename, unnamed using LinearAlgebra: norm, normalize! -using MatrixAlgebraKit: project_hermitian, qr_compact, svd_trunc +using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc using NamedGraphs: boundary_edges -using TensorAlgebra.MatrixAlgebra: sqrth_invsqrth_safe, sqrth_safe +using TensorAlgebra.MatrixAlgebra: invsqrth_safe, sqrth_safe +using TensorAlgebra: isdual, matricize, twist!, unmatricize + +# Asymmetric (Gram) root of a Hermitian positive semidefinite matrix, as the pair +# `(root, inv_root)`: `root' * root == m`, and `inv_root * root` is the identity on `m`'s +# support. `eigh_full` is a square decomposition, so both factors are square. +function gram_root_invroot(m::AbstractMatrix; kwargs...) + d, u = eigh_full(m) + return sqrth_safe(d; kwargs...) * u', u * invsqrth_safe(d; kwargs...) +end + +# The same root for a named tensor split into `outnames` and `innames`, as the pair `(x, y)`: +# `x * ψ` absorbs it into a state tensor `ψ` and `y * ·` un-absorbs it. +# +# `y` inverts `x` as a matrix, and `contract` is that matrix product plus a twist, so the two +# have different identity elements. Twisting `y` on its dual bond axes, the same axes `contract` +# twists, makes it the inverse under `contract`. `twist!` is the identity on storage carrying no +# sector data, so this is a no-op outside the fermionic case. +# +# The message is matricized ket to bra, which keeps the absorbed wavefunction ket-like, and the +# eigendecomposition hands back the rank space already carrying that orientation. +function gram_root_invroot(t::AbstractNamedTensor, outnames, innames; kwargs...) + bondinds = Tuple(conj.(to_inds(t, innames))) + root, inv_root = gram_root_invroot(matricize(t, outnames, innames); kwargs...) + rankind = Index(axes(root, 1)) + y = unmatricize(inv_root, bondinds, (rankind,)) + twist!(y, filter(isdual, bondinds)) + return unmatricize(root, (rankind,), bondinds), y +end + +function gram_root_invroot(m::NamedTensorOperator; kwargs...) + return gram_root_invroot(state(m), outputnames(m), inputnames(m); kwargs...) +end # === Top-level user entry point === @@ -252,14 +286,12 @@ function apply_gate_bp_nsite!( v1, v2 = vertices edges_in = boundary_edges(state, vertices; dir = :in) roots_v1 = - [sqrth_invsqrth_safe(project_hermitian(env[e])) for e in edges_in if dst(e) == v1] + [gram_root_invroot(project_hermitian(env[e])) for e in edges_in if dst(e) == v1] roots_v2 = - [sqrth_invsqrth_safe(project_hermitian(env[e])) for e in edges_in if dst(e) == v2] - sqrt_messages_v1, invsqrt_messages_v1 = first.(roots_v1), last.(roots_v1) - sqrt_messages_v2, invsqrt_messages_v2 = first.(roots_v2), last.(roots_v2) + [gram_root_invroot(project_hermitian(env[e])) for e in edges_in if dst(e) == v2] - ψ_v1 = foldl((ψ, m) -> apply(m, ψ), sqrt_messages_v1; init = state[v1]) - ψ_v2 = foldl((ψ, m) -> apply(m, ψ), sqrt_messages_v2; init = state[v2]) + ψ_v1 = foldl((ψ, (x, _)) -> x * ψ, roots_v1; init = state[v1]) + ψ_v2 = foldl((ψ, (x, _)) -> x * ψ, roots_v2; init = state[v2]) Q_v1, R_v1 = qr_compact(ψ_v1, setdiff(names(ψ_v1), names(ψ_v2), names(op))) Q_v2, R_v2 = qr_compact(ψ_v2, setdiff(names(ψ_v2), names(ψ_v1), names(op))) @@ -273,8 +305,8 @@ function apply_gate_bp_nsite!( R_v1 = rename(U_v1 * sqrt_S, name_v2 => name_v1) R_v2 = sqrt_S * U_v2 - dest[v1] = foldl((ψ, m) -> apply(m, ψ), invsqrt_messages_v1; init = Q_v1 * R_v1) - dest[v2] = foldl((ψ, m) -> apply(m, ψ), invsqrt_messages_v2; init = Q_v2 * R_v2) + dest[v1] = foldl((ψ, (_, y)) -> y * ψ, roots_v1; init = Q_v1 * R_v1) + dest[v2] = foldl((ψ, (_, y)) -> y * ψ, roots_v2; init = Q_v2 * R_v2) env[v1 => v2] = operator( rename(conj(R_v1), name_v1 => name_v2) * R_v1, (name_v2,), (name_v1,) diff --git a/test/Project.toml b/test/Project.toml index d761064c..9abc495b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -31,7 +31,7 @@ DataGraphs = "0.6" Dictionaries = "0.4.5" GradedArrays = "0.17" Graphs = "1.13.1" -ITensorBase = "0.16" +ITensorBase = "0.16.1" ITensorNetworksNext = "0.11" ITensorPkgSkeleton = "0.3.42" MatrixAlgebraKit = "0.6" @@ -40,6 +40,6 @@ OMEinsumContractionOrders = "1" QuadGK = "2.11.2" Random = "1.10" StableRNGs = "1" -TensorAlgebra = "0.23" +TensorAlgebra = "0.23.1" TensorKitSectors = "0.3" Test = "1.10" From 81f9b342c804ba856c8c9e97bc025e53149ea5dc Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 7 Oct 2026 09:15:03 -0400 Subject: [PATCH 2/5] Name the bond-message gauge, and project it Hermitian in one place Calls the operator-level root `message_gauge`, since that is the layer that knows its argument is a bond message, and moves the Hermitian projection the call sites were each doing into it. The two generic layers keep the `gram_root_invroot` name, which describes the factorization rather than the use. Co-Authored-By: Claude Opus 5 (1M context) --- src/apply/apply_operators.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl index 5eaf68b9..9bee0a46 100644 --- a/src/apply/apply_operators.jl +++ b/src/apply/apply_operators.jl @@ -38,8 +38,12 @@ function gram_root_invroot(t::AbstractNamedTensor, outnames, innames; kwargs...) return unmatricize(root, (rankind,), bondinds), y end -function gram_root_invroot(m::NamedTensorOperator; kwargs...) - return gram_root_invroot(state(m), outputnames(m), inputnames(m); kwargs...) +# The gauge a bond message induces: its output names are the bra side and its input names the +# ket side, which is the split the root is taken over. A message is Hermitian only up to +# numerical noise, and the root needs it exactly so. +function message_gauge(m::NamedTensorOperator; kwargs...) + h = project_hermitian(m) + return gram_root_invroot(state(h), outputnames(h), inputnames(h); kwargs...) end # === Top-level user entry point === @@ -286,9 +290,9 @@ function apply_gate_bp_nsite!( v1, v2 = vertices edges_in = boundary_edges(state, vertices; dir = :in) roots_v1 = - [gram_root_invroot(project_hermitian(env[e])) for e in edges_in if dst(e) == v1] + [message_gauge(env[e]) for e in edges_in if dst(e) == v1] roots_v2 = - [gram_root_invroot(project_hermitian(env[e])) for e in edges_in if dst(e) == v2] + [message_gauge(env[e]) for e in edges_in if dst(e) == v2] ψ_v1 = foldl((ψ, (x, _)) -> x * ψ, roots_v1; init = state[v1]) ψ_v2 = foldl((ψ, (x, _)) -> x * ψ, roots_v2; init = state[v2]) From 87abdc02b8312978b187222598ff79bb90747c90 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 7 Oct 2026 12:24:58 -0400 Subject: [PATCH 3/5] Resolve ITensorBase from the registry now that 0.16.1 is released Co-Authored-By: Claude Opus 5 (1M context) --- Project.toml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Project.toml b/Project.toml index deff89b0..c614f208 100644 --- a/Project.toml +++ b/Project.toml @@ -27,10 +27,6 @@ TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" [weakdeps] OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715" -[sources.ITensorBase] -rev = "mf/named-twist" -url = "https://github.com/ITensor/ITensorBase.jl" - [extensions] ITensorNetworksNextOMEinsumContractionOrdersExt = "OMEinsumContractionOrders" From 68998db0173b6bf0315958370aebb23cceb3e86a Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 7 Oct 2026 12:49:22 -0400 Subject: [PATCH 4/5] Call the whole gauge family message_gauge Keeps one name for the operation across its matrix, named-tensor and operator layers, rather than naming the lower two after the factorization that implements them. Co-Authored-By: Claude Opus 5 (1M context) --- src/apply/apply_operators.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl index 9bee0a46..13db7916 100644 --- a/src/apply/apply_operators.jl +++ b/src/apply/apply_operators.jl @@ -14,7 +14,7 @@ using TensorAlgebra: isdual, matricize, twist!, unmatricize # Asymmetric (Gram) root of a Hermitian positive semidefinite matrix, as the pair # `(root, inv_root)`: `root' * root == m`, and `inv_root * root` is the identity on `m`'s # support. `eigh_full` is a square decomposition, so both factors are square. -function gram_root_invroot(m::AbstractMatrix; kwargs...) +function message_gauge(m::AbstractMatrix; kwargs...) d, u = eigh_full(m) return sqrth_safe(d; kwargs...) * u', u * invsqrth_safe(d; kwargs...) end @@ -29,21 +29,21 @@ end # # The message is matricized ket to bra, which keeps the absorbed wavefunction ket-like, and the # eigendecomposition hands back the rank space already carrying that orientation. -function gram_root_invroot(t::AbstractNamedTensor, outnames, innames; kwargs...) +function message_gauge(t::AbstractNamedTensor, outnames, innames; kwargs...) bondinds = Tuple(conj.(to_inds(t, innames))) - root, inv_root = gram_root_invroot(matricize(t, outnames, innames); kwargs...) + root, inv_root = message_gauge(matricize(t, outnames, innames); kwargs...) rankind = Index(axes(root, 1)) y = unmatricize(inv_root, bondinds, (rankind,)) twist!(y, filter(isdual, bondinds)) return unmatricize(root, (rankind,), bondinds), y end -# The gauge a bond message induces: its output names are the bra side and its input names the -# ket side, which is the split the root is taken over. A message is Hermitian only up to -# numerical noise, and the root needs it exactly so. +# A message's output names are the bra side and its input names the ket side, which is the split +# the gauge is taken over. A message is Hermitian only up to numerical noise, and the root needs +# it exactly so. function message_gauge(m::NamedTensorOperator; kwargs...) h = project_hermitian(m) - return gram_root_invroot(state(h), outputnames(h), inputnames(h); kwargs...) + return message_gauge(state(h), outputnames(h), inputnames(h); kwargs...) end # === Top-level user entry point === From 24c7419fe9ec65037eb7c157da7181fbab8c996f Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 7 Oct 2026 13:18:31 -0400 Subject: [PATCH 5/5] Pass the gauge its indices rather than its names Takes the bond indices directly, through `inputinds`, instead of looking names up and dualizing them by hand at the one call site that knew the convention. Co-Authored-By: Claude Opus 5 (1M context) --- src/apply/apply_operators.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/apply/apply_operators.jl b/src/apply/apply_operators.jl index 13db7916..6736d655 100644 --- a/src/apply/apply_operators.jl +++ b/src/apply/apply_operators.jl @@ -3,8 +3,8 @@ using AlgorithmsInterface: AlgorithmsInterface as AI using Base: @kwdef using Graphs: dst, src, vertices using ITensorBase: AbstractITensor, AbstractNamedTensor, ITensor, Index, NamedTensor, - NamedTensorOperator, apply, inputinds, inputnames, name, names, operator, outputnames, - rename, sim, state, to_inds, uniquename, unnamed + NamedTensorOperator, apply, inputinds, inputnames, name, names, operator, outputinds, + rename, sim, state, uniquename, unnamed using LinearAlgebra: norm, normalize! using MatrixAlgebraKit: eigh_full, project_hermitian, qr_compact, svd_trunc using NamedGraphs: boundary_edges @@ -19,7 +19,7 @@ function message_gauge(m::AbstractMatrix; kwargs...) return sqrth_safe(d; kwargs...) * u', u * invsqrth_safe(d; kwargs...) end -# The same root for a named tensor split into `outnames` and `innames`, as the pair `(x, y)`: +# The same root for a named tensor split into `outinds` and `bondinds`, as the pair `(x, y)`: # `x * ψ` absorbs it into a state tensor `ψ` and `y * ·` un-absorbs it. # # `y` inverts `x` as a matrix, and `contract` is that matrix product plus a twist, so the two @@ -29,9 +29,8 @@ end # # The message is matricized ket to bra, which keeps the absorbed wavefunction ket-like, and the # eigendecomposition hands back the rank space already carrying that orientation. -function message_gauge(t::AbstractNamedTensor, outnames, innames; kwargs...) - bondinds = Tuple(conj.(to_inds(t, innames))) - root, inv_root = message_gauge(matricize(t, outnames, innames); kwargs...) +function message_gauge(t::AbstractNamedTensor, outinds, bondinds; kwargs...) + root, inv_root = message_gauge(matricize(t, outinds, bondinds); kwargs...) rankind = Index(axes(root, 1)) y = unmatricize(inv_root, bondinds, (rankind,)) twist!(y, filter(isdual, bondinds)) @@ -43,7 +42,7 @@ end # it exactly so. function message_gauge(m::NamedTensorOperator; kwargs...) h = project_hermitian(m) - return message_gauge(state(h), outputnames(h), inputnames(h); kwargs...) + return message_gauge(state(h), Tuple(outputinds(h)), Tuple(inputinds(h)); kwargs...) end # === Top-level user entry point ===