From 3dd60136f8ae9002fadd43f1e432e5143113acfa Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 20 Jul 2026 12:01:34 -0400 Subject: [PATCH 1/3] Add operator forms of the Hermitian square-root factorizations Adds operator (`NamedTensorOperator`) forms of `project_hermitian`, `sqrth_safe`, `invsqrth_safe`, and `sqrth_invsqrth_safe`, completing the operator-form factorization family alongside the existing `gram_eigh_full` operator forms. Each takes its codomain/domain bipartition from the operator, so callers can compose it with `transpose` to project Hermitian in one bipartition and take the root in the transposed one. `project_hermitian` returns an operator and the roots return bare named arrays, matching the tensor-level forms. --- Project.toml | 2 +- src/namedtensoroperator.jl | 16 ++++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0af58ad..17f558c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorBase" uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" -version = "0.12.2" +version = "0.12.3" authors = ["ITensor developers and contributors"] [workspace] diff --git a/src/namedtensoroperator.jl b/src/namedtensoroperator.jl index f5c3526..41356bf 100644 --- a/src/namedtensoroperator.jl +++ b/src/namedtensoroperator.jl @@ -438,6 +438,22 @@ for f in (:gram_eigh_full, :gram_eigh_full_with_pinv) end end +# Operator forms of the Hermitian-square-root family: the codomain/domain partition is +# taken from the operator, so callers can compose with `transpose` to choose the +# bipartition — e.g. project Hermitian in one bipartition and take the square root in the +# transposed one, which induces the fermionic braid sign on the odd-parity sector. +# `project_hermitian` returns an operator (its Hermitian part is again a bond operator); the +# roots return bare named arrays like `gram_eigh_full`, being terminal factorization outputs. +function MAK.project_hermitian(a::NamedTensorOperator; kwargs...) + h = MAK.project_hermitian(state(a), codomainnames(a), domainnames(a); kwargs...) + return operator(h, codomainnames(a), domainnames(a)) +end +for f in (:sqrth_safe, :invsqrth_safe, :sqrth_invsqrth_safe) + @eval function MA.$f(a::NamedTensorOperator; kwargs...) + return MA.$f(state(a), codomainnames(a), domainnames(a); kwargs...) + end +end + """ Base.one(op::NamedTensorOperator) -> Id From 10d02b1b4eed04991ca9cc3e9131fc12ed950d21 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 20 Jul 2026 12:11:32 -0400 Subject: [PATCH 2/3] Add tests for the operator-form Hermitian square roots Covers project_hermitian, sqrth_safe, invsqrth_safe, and sqrth_invsqrth_safe on NamedTensorOperator, matching the gram_eigh_full test pattern. --- src/namedtensoroperator.jl | 10 ++++------ test/test_operator.jl | 34 +++++++++++++++++++++++++++++++++- 2 files changed, 37 insertions(+), 7 deletions(-) diff --git a/src/namedtensoroperator.jl b/src/namedtensoroperator.jl index 41356bf..9aa6501 100644 --- a/src/namedtensoroperator.jl +++ b/src/namedtensoroperator.jl @@ -438,12 +438,10 @@ for f in (:gram_eigh_full, :gram_eigh_full_with_pinv) end end -# Operator forms of the Hermitian-square-root family: the codomain/domain partition is -# taken from the operator, so callers can compose with `transpose` to choose the -# bipartition — e.g. project Hermitian in one bipartition and take the square root in the -# transposed one, which induces the fermionic braid sign on the odd-parity sector. -# `project_hermitian` returns an operator (its Hermitian part is again a bond operator); the -# roots return bare named arrays like `gram_eigh_full`, being terminal factorization outputs. +# The codomain/domain bipartition is taken from the operator, so callers can compose with +# `transpose` to choose it. `project_hermitian` returns an operator (its Hermitian part is +# again a bond operator); the roots return bare named arrays like `gram_eigh_full`, being +# terminal factorization outputs. function MAK.project_hermitian(a::NamedTensorOperator; kwargs...) h = MAK.project_hermitian(state(a), codomainnames(a), domainnames(a); kwargs...) return operator(h, codomainnames(a), domainnames(a)) diff --git a/test/test_operator.jl b/test/test_operator.jl index cfb08a7..a9bf18b 100644 --- a/test/test_operator.jl +++ b/test/test_operator.jl @@ -2,9 +2,11 @@ using ITensorBase: ITensorBase as NDA, NamedTensor, NamedTensorOperator, apply, codomainnames, dimnames, domainnames, id, nameddims, namedoneto, operator, product, replacedimnames, similar_operator, state, unname, unnamed using LinearAlgebra: I, norm +using MatrixAlgebraKit: project_hermitian using Random: Random using StableRNGs: StableRNG -using TensorAlgebra.MatrixAlgebra: gram_eigh_full, gram_eigh_full_with_pinv +using TensorAlgebra.MatrixAlgebra: + gram_eigh_full, gram_eigh_full_with_pinv, invsqrth_safe, sqrth_invsqrth_safe, sqrth_safe using TensorAlgebra: matricize using Test: @test, @test_throws, @testset @@ -229,3 +231,33 @@ end @test Xp2 * Xp2' ≈ A @test Yp2 * Xp2 ≈ I(n) end + +@testset "Hermitian square roots on NamedTensorOperator" begin + n = 5 + B = randn(n, n) + A = B * B' # Hermitian PSD + M_nda = nameddims(A, ("ket", "bra")) + M_op = operator(M_nda, ["ket"], ["bra"]) + + # `project_hermitian` keeps the operator structure; a non-Hermitian input maps to its + # Hermitian part. + H_op = project_hermitian(M_op) + @test H_op isa NamedTensorOperator + @test codomainnames(H_op) == codomainnames(M_op) + @test domainnames(H_op) == domainnames(M_op) + @test H_op ≈ M_op + N_op = operator(nameddims(B, ("ket", "bra")), ["ket"], ["bra"]) + @test project_hermitian(N_op) ≈ + operator(nameddims((B + B') / 2, ("ket", "bra")), ["ket"], ["bra"]) + + # The roots forward to the named-array entry, so match its data and shape. + @test parent(sqrth_safe(M_op)) ≈ parent(sqrth_safe(M_nda, ("ket",), ("bra",))) + + P = parent(sqrth_safe(M_op)) + @test P * P' ≈ A + @test parent(invsqrth_safe(M_op)) * P ≈ I(n) + + P2, Pinv2 = sqrth_invsqrth_safe(M_op) + @test parent(P2) * parent(P2)' ≈ A + @test parent(P2) * parent(Pinv2) ≈ I(n) +end From 8e904231e318196b8321ca208fec248b4c4c34f2 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 20 Jul 2026 12:43:06 -0400 Subject: [PATCH 3/3] Return operators from the operator-form square roots A square root of a bond operator is again a bond operator, so sqrth_safe, invsqrth_safe, and sqrth_invsqrth_safe now return NamedTensorOperator with the same codomain/domain as the input rather than bare named arrays. --- src/namedtensoroperator.jl | 14 ++++++++------ test/test_operator.jl | 34 ++++++++++++++++++---------------- 2 files changed, 26 insertions(+), 22 deletions(-) diff --git a/src/namedtensoroperator.jl b/src/namedtensoroperator.jl index 9aa6501..77cc09b 100644 --- a/src/namedtensoroperator.jl +++ b/src/namedtensoroperator.jl @@ -438,19 +438,21 @@ for f in (:gram_eigh_full, :gram_eigh_full_with_pinv) end end -# The codomain/domain bipartition is taken from the operator, so callers can compose with -# `transpose` to choose it. `project_hermitian` returns an operator (its Hermitian part is -# again a bond operator); the roots return bare named arrays like `gram_eigh_full`, being -# terminal factorization outputs. function MAK.project_hermitian(a::NamedTensorOperator; kwargs...) h = MAK.project_hermitian(state(a), codomainnames(a), domainnames(a); kwargs...) return operator(h, codomainnames(a), domainnames(a)) end -for f in (:sqrth_safe, :invsqrth_safe, :sqrth_invsqrth_safe) +for f in (:sqrth_safe, :invsqrth_safe) @eval function MA.$f(a::NamedTensorOperator; kwargs...) - return MA.$f(state(a), codomainnames(a), domainnames(a); kwargs...) + x = MA.$f(state(a), codomainnames(a), domainnames(a); kwargs...) + return operator(x, codomainnames(a), domainnames(a)) end end +function MA.sqrth_invsqrth_safe(a::NamedTensorOperator; kwargs...) + x, y = MA.sqrth_invsqrth_safe(state(a), codomainnames(a), domainnames(a); kwargs...) + return operator(x, codomainnames(a), domainnames(a)), + operator(y, codomainnames(a), domainnames(a)) +end """ Base.one(op::NamedTensorOperator) -> Id diff --git a/test/test_operator.jl b/test/test_operator.jl index a9bf18b..b3d3bba 100644 --- a/test/test_operator.jl +++ b/test/test_operator.jl @@ -214,11 +214,10 @@ end n = 5 B = randn(n, n) A = B * B' # Hermitian PSD - M_nda = nameddims(A, ("ket", "bra")) - M_op = operator(M_nda, ["ket"], ["bra"]) + M_op = operator(A, ["ket"], ["bra"]) X_op = gram_eigh_full(M_op) - X_arr = gram_eigh_full(M_nda, ("ket",), ("bra",)) + X_arr = gram_eigh_full(nameddims(A, ("ket", "bra")), ("ket",), ("bra",)) # Operator entry forwards to the named-array entry: same data, same shape. @test size(parent(X_op)) == size(parent(X_arr)) @@ -236,8 +235,7 @@ end n = 5 B = randn(n, n) A = B * B' # Hermitian PSD - M_nda = nameddims(A, ("ket", "bra")) - M_op = operator(M_nda, ["ket"], ["bra"]) + M_op = operator(A, ["ket"], ["bra"]) # `project_hermitian` keeps the operator structure; a non-Hermitian input maps to its # Hermitian part. @@ -246,18 +244,22 @@ end @test codomainnames(H_op) == codomainnames(M_op) @test domainnames(H_op) == domainnames(M_op) @test H_op ≈ M_op - N_op = operator(nameddims(B, ("ket", "bra")), ["ket"], ["bra"]) - @test project_hermitian(N_op) ≈ - operator(nameddims((B + B') / 2, ("ket", "bra")), ["ket"], ["bra"]) - - # The roots forward to the named-array entry, so match its data and shape. - @test parent(sqrth_safe(M_op)) ≈ parent(sqrth_safe(M_nda, ("ket",), ("bra",))) + @test project_hermitian(operator(B, ["ket"], ["bra"])) ≈ + operator((B + B') / 2, ["ket"], ["bra"]) + + # The roots are again bond operators, with the same codomain/domain as the input. + for X in (sqrth_safe(M_op), invsqrth_safe(M_op), sqrth_invsqrth_safe(M_op)...) + @test X isa NamedTensorOperator + @test codomainnames(X) == codomainnames(M_op) + @test domainnames(X) == domainnames(M_op) + end - P = parent(sqrth_safe(M_op)) + P = unnamed(state(sqrth_safe(M_op))) @test P * P' ≈ A - @test parent(invsqrth_safe(M_op)) * P ≈ I(n) + @test unnamed(state(invsqrth_safe(M_op))) * P ≈ I(n) - P2, Pinv2 = sqrth_invsqrth_safe(M_op) - @test parent(P2) * parent(P2)' ≈ A - @test parent(P2) * parent(Pinv2) ≈ I(n) + Psqrt, Pinv = sqrth_invsqrth_safe(M_op) + Pmat = unnamed(state(Psqrt)) + @test Pmat * Pmat' ≈ A + @test Pmat * unnamed(state(Pinv)) ≈ I(n) end