Add operator forms of the Hermitian square-root factorizations

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 <support@itensor.org> and contributors"]
 
 [workspace]

src/namedtensoroperator.jl

@@ -438,6 +438,22 @@ for f in (:gram_eigh_full, :gram_eigh_full_with_pinv)
     end
 end
 
+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)
+    @eval function MA.$f(a::NamedTensorOperator; 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
 

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
 
@@ -212,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))
 
@@ -229,3 +230,36 @@ 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_op = operator(A, ["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
+    @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 = unnamed(state(sqrth_safe(M_op)))
+    @test P * P' ≈ A
+    @test unnamed(state(invsqrth_safe(M_op))) * P ≈ 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