Follow TensorAlgebra v0.21

Project.toml

@@ -1,6 +1,6 @@
 name = "ITensorBase"
 uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
-version = "0.14.1"
+version = "0.14.2"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
 
 [workspace]
@@ -45,14 +45,14 @@ Adapt = "4.1.1"
 ArrayLayouts = "1.11"
 Combinatorics = "1"
 ConstructionBase = "1.6"
-GradedArrays = "0.16"
+GradedArrays = "0.16.4"
 LinearAlgebra = "1.10"
 MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5, 0.6"
 Mooncake = "0.4.202, 0.5"
 OMEinsumContractionOrders = "1.3"
 Random = "1.10"
 SimpleTraits = "0.9.4"
-TensorAlgebra = "0.20"
+TensorAlgebra = "0.21"
 TensorKit = "0.17"
 TensorKitSectors = "0.3.9"
 TermInterface = "2"

docs/Project.toml

@@ -16,5 +16,5 @@ ITensorBase = "0.14"
 ITensorFormatter = "0.2.27"
 Literate = "2"
 MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5, 0.6"
-TensorAlgebra = "0.20"
+TensorAlgebra = "0.21"
 Test = "1.10"

src/namedtensoroperator.jl

@@ -517,90 +517,9 @@ for f in MATRIX_FUNCTIONS
     end
 end
 
-# Operator entries for the gram factorizations defined in `tensoralgebra.jl`.
+# Operator entries for the Hermitian factorizations defined in `tensoralgebra.jl`.
 # Placed here because `NamedTensorOperator` is defined in this file, which comes
 # after `tensoralgebra.jl` in the include order.
-#
-# Per-method docstrings are factored out into `const` strings and attached
-# inside the `@eval` loop via `@doc`. This keeps the loop body uniform when
-# methods need distinct user-facing docs (including jldoctest examples) that
-# don't share enough structure to warrant `$($f)`-interpolation.
-
-const _gram_eigh_full_operator_docstring = """
-    TensorAlgebra.MatrixAlgebra.gram_eigh_full(a::NamedTensorOperator; kwargs...) -> x
-
-Gram factorization of a Hermitian positive semi-definite named operator
-`a`, returning `x` such that `x * x_cod ≈ state(a)`, where `x_cod` is
-`conj(x)` with its input dimension names replaced by the corresponding
-output names of `a`. `x` carries `a`'s input dimension names and a
-fresh trailing rank name. The output and input partition is taken from
-`outputnames(a)` and `inputnames(a)`.
-
-`kwargs` are forwarded to `TensorAlgebra.MatrixAlgebra.gram_eigh_full` on the
-underlying named array (e.g. `atol`, `rtol`).
-
-# Examples
-
-```jldoctest
-julia> using ITensorBase: namedoneto, operator, replacedimnames, state
-
-julia> using TensorAlgebra.MatrixAlgebra: gram_eigh_full
-
-julia> i, j, k, l, aux = namedoneto.((2, 2, 2, 2, 8), ("i", "j", "k", "l", "aux"));
-
-julia> b = randn(aux, i, k);
-
-julia> a = operator(conj(b) * replacedimnames(b, "i" => "j", "k" => "l"), ("i", "k"), ("j", "l"));
-
-julia> x = gram_eigh_full(a);
-
-julia> replacedimnames(x, "j" => "i", "l" => "k") * conj(x) ≈ state(a)
-true
-```
-"""
-
-const _gram_eigh_full_with_pinv_operator_docstring = """
-    TensorAlgebra.MatrixAlgebra.gram_eigh_full_with_pinv(a::NamedTensorOperator; kwargs...) -> x, y
-
-Like `TensorAlgebra.MatrixAlgebra.gram_eigh_full`, but additionally returns a
-named array `y` that is a left inverse of `x`: `y * x ≈ I` on the
-rank subspace (equal to the identity when `a` is full rank). The
-output and input partition is taken from `outputnames(a)` and
-`inputnames(a)`.
-
-# Examples
-
-```jldoctest
-julia> using LinearAlgebra: I
-
-julia> using ITensorBase: unname, dimnames, namedoneto, operator, replacedimnames
-
-julia> using TensorAlgebra.MatrixAlgebra: gram_eigh_full_with_pinv
-
-julia> i, j, k, l, aux = namedoneto.((2, 2, 2, 2, 8), ("i", "j", "k", "l", "aux"));
-
-julia> b = randn(aux, i, k);
-
-julia> a = operator(conj(b) * replacedimnames(b, "i" => "j", "k" => "l"), ("i", "k"), ("j", "l"));
-
-julia> x, y = gram_eigh_full_with_pinv(a);
-
-julia> rname = only(setdiff(dimnames(x), ("j", "l")));
-
-julia> reshape(unname(y, (rname, "j", "l")), :, 4) *
-       reshape(unname(x, ("j", "l", rname)), 4, :) ≈ I
-true
-```
-"""
-
-for f in (:gram_eigh_full, :gram_eigh_full_with_pinv)
-    doc_sym = Symbol("_", f, "_operator_docstring")
-    @eval begin
-        @doc $doc_sym function MA.$f(a::NamedTensorOperator; kwargs...)
-            return MA.$f(state(a), outputnames(a), inputnames(a); kwargs...)
-        end
-    end
-end
 
 function MAK.project_hermitian(a::NamedTensorOperator; kwargs...)
     h = MAK.project_hermitian(state(a), outputnames(a), inputnames(a); kwargs...)

src/tensoralgebra.jl

@@ -120,7 +120,7 @@ end
 function matricize_nameddims(na::AbstractNamedTensor, fusions::Vararg{Pair, 2})
     group1, group2 = first.(fusions)
     perm_codomain, perm_domain = nameperm(na, group1, group2)
-    a_fused = TA.matricizeperm(unnamed(na), perm_codomain, perm_domain)
+    a_fused = TA.matricize(unnamed(na), perm_codomain, perm_domain)
     return nameddims(a_fused, last.(fusions))
 end
 
@@ -467,111 +467,6 @@ function right_null_nameddims(a::AbstractNamedTensor, dimnames_codomain; kwargs.
     return MAK.right_null(a, codomain, domain; kwargs...)
 end
 
-"""
-    TensorAlgebra.MatrixAlgebra.gram_eigh_full(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> x
-
-Gram factorization of a Hermitian positive semi-definite named array `a`,
-returning `x` such that `a ≈ x * x_cod`, where `x_cod` is `conj(x)` with
-its domain dimension names replaced by the corresponding codomain names.
-`x` carries the domain dimension names of `a` (matching the convention
-that the stored factor labels a vector in `a`'s input space) and a fresh
-trailing rank name.
-
-`kwargs` are forwarded to `TensorAlgebra.gram_eigh_full` on the underlying
-unnamed array (e.g. `atol`, `rtol`).
-
-# Examples
-
-```jldoctest
-julia> using ITensorBase: dimnames, namedoneto, replacedimnames
-
-julia> using TensorAlgebra.MatrixAlgebra: gram_eigh_full
-
-julia> i, j, k, l, aux = namedoneto.((2, 2, 2, 2, 8), ("i", "j", "k", "l", "aux"));
-
-julia> b = randn(aux, i, k);
-
-julia> a = conj(b) * replacedimnames(b, "i" => "j", "k" => "l");
-
-julia> x = gram_eigh_full(a, (i, k), (j, l));
-
-julia> replacedimnames(x, "j" => "i", "l" => "k") * conj(x) ≈ a
-true
-```
-"""
-function MA.gram_eigh_full(
-        a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...
-    )
-    return gram_eigh_full_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...)
-end
-function gram_eigh_full_nameddims(
-        a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...
-    )
-    codomain = name.(dimnames_codomain)
-    domain = name.(dimnames_domain)
-    x_unnamed = TA.gram_eigh_full(unnamed(a), dimnames(a), codomain, domain; kwargs...)
-    name_x = uniquename(dimnametype(a))
-    dimnames_x = (domain..., name_x)
-    return nameddims(x_unnamed, dimnames_x)
-end
-
-"""
-    TensorAlgebra.MatrixAlgebra.gram_eigh_full_with_pinv(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> x, y
-
-Like `TensorAlgebra.MatrixAlgebra.gram_eigh_full`, but additionally returns a
-named array `y` that is a left inverse of `x`: `y * x ≈ I` on the rank
-subspace (equal to the identity when `a` is full rank). `x` has the
-rank-name last, `y` has it first, both sharing the domain dimension
-names of `a`.
-
-# Examples
-
-```jldoctest
-julia> using LinearAlgebra: I
-
-julia> using ITensorBase: unname, dimnames, namedoneto, replacedimnames
-
-julia> using TensorAlgebra.MatrixAlgebra: gram_eigh_full_with_pinv
-
-julia> i, j, k, l, aux = namedoneto.((2, 2, 2, 2, 8), ("i", "j", "k", "l", "aux"));
-
-julia> b = randn(aux, i, k);
-
-julia> a = conj(b) * replacedimnames(b, "i" => "j", "k" => "l");
-
-julia> x, y = gram_eigh_full_with_pinv(a, (i, k), (j, l));
-
-julia> replacedimnames(x, "j" => "i", "l" => "k") * conj(x) ≈ a
-true
-
-julia> rname = only(setdiff(dimnames(x), ("j", "l")));
-
-julia> reshape(unname(y, (rname, "j", "l")), :, 4) *
-       reshape(unname(x, ("j", "l", rname)), 4, :) ≈ I
-true
-```
-"""
-function MA.gram_eigh_full_with_pinv(
-        a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...
-    )
-    return gram_eigh_full_with_pinv_nameddims(
-        a, dimnames_codomain, dimnames_domain; kwargs...
-    )
-end
-function gram_eigh_full_with_pinv_nameddims(
-        a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...
-    )
-    codomain = name.(dimnames_codomain)
-    domain = name.(dimnames_domain)
-    x_unnamed, y_unnamed = TA.gram_eigh_full_with_pinv(
-        unnamed(a), dimnames(a), codomain, domain; kwargs...
-    )
-    name_xy = uniquename(dimnametype(a))
-    dimnames_x = (domain..., name_xy)
-    dimnames_y = (name_xy, domain...)
-    return nameddims(x_unnamed, dimnames_x), nameddims(y_unnamed, dimnames_y)
-end
-
 """
     TensorAlgebra.MatrixAlgebra.sqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p
 

test/Project.toml

@@ -32,7 +32,7 @@ AbstractTrees = "0.4.5"
 Adapt = "4"
 Aqua = "0.8.9"
 Combinatorics = "1"
-GradedArrays = "0.16"
+GradedArrays = "0.16.4"
 ITensorBase = "0.14"
 ITensorPkgSkeleton = "0.3.42"
 JLArrays = "0.2, 0.3"
@@ -44,7 +44,7 @@ Random = "1.10"
 SafeTestsets = "0.1"
 StableRNGs = "1"
 Suppressor = "0.2"
-TensorAlgebra = "0.20"
+TensorAlgebra = "0.21"
 TensorKit = "0.17"
 TensorKitSectors = "0.3.9"
 TermInterface = "2"

test/test_operator.jl

@@ -7,8 +7,7 @@ using LinearAlgebra: I, norm
 using MatrixAlgebraKit: project_hermitian
 using Random: Random, randn
 using StableRNGs: StableRNG
-using TensorAlgebra.MatrixAlgebra:
-    gram_eigh_full, gram_eigh_full_with_pinv, invsqrth_safe, sqrth_invsqrth_safe, sqrth_safe
+using TensorAlgebra.MatrixAlgebra: invsqrth_safe, sqrth_invsqrth_safe, sqrth_safe
 using TensorAlgebra: matricize
 using Test: @test, @test_throws, @testset
 
@@ -420,27 +419,6 @@ end
     @test issetequal(dimnames(Av), ("a'",))
 end
 
-@testset "gram_eigh_full on NamedTensorOperator" begin
-    n = 5
-    B = randn(n, n)
-    A = B * B'  # Hermitian PSD
-    M_op = operator(A, ["ket"], ["bra"])
-
-    X_op = gram_eigh_full(M_op)
-    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))
-
-    Xp = parent(X_op)
-    @test Xp * Xp' ≈ A
-
-    X2, Y2 = gram_eigh_full_with_pinv(M_op)
-    Xp2 = parent(X2)
-    Yp2 = parent(Y2)
-    @test Xp2 * Xp2' ≈ A
-    @test Yp2 * Xp2 ≈ I(n)
-end
-
 @testset "Hermitian square roots on NamedTensorOperator" begin
     n = 5
     B = randn(n, n)

test/test_tensoralgebra.jl

@@ -1,10 +1,9 @@
 using ITensorBase: ITensorBase, Index, dimnames, id, inds, name, namedoneto, operator,
-    prime, replacedimnames, uniquename, unname, unnamed
-using LinearAlgebra: LinearAlgebra, norm, tr
+    prime, replacedimnames, unname, unnamed
+using LinearAlgebra: norm, tr
 using MatrixAlgebraKit: left_null, left_orth, left_polar, lq_compact, lq_full, qr_compact,
     qr_full, right_null, right_orth, right_polar, svd_compact, svd_trunc, svd_vals
 using StableRNGs: StableRNG
-using TensorAlgebra.MatrixAlgebra: gram_eigh_full, gram_eigh_full_with_pinv
 using TensorAlgebra: TensorAlgebra, contract, directsum, matricize, project, trivialrange,
     unchecked_project, unmatricize
 using Test: @test, @test_broken, @testset
@@ -172,36 +171,6 @@ using Test: @test, @test_broken, @testset
             @test norm(n * a) ≈ 0
         end
     end
-    @testset "gram_eigh_full" begin
-        # Build a Hermitian PSD a ≈ conj(b) * b over an aux dim, with codomain
-        # (i, k) and domain (j, l) sharing the same axis lengths.
-        i, j, k, l, aux = namedoneto.((2, 2, 2, 2, 5), ("i", "j", "k", "l", "aux"))
-        b = randn(elt, aux, i, k)
-        # conj(b) * b with the non-conjugated copy's (i, k) relabeled to
-        # (j, l) to form the operator-shaped Hermitian a ≈ X * X'.
-        b_dom = replacedimnames(b, "i" => "j", "k" => "l")
-        a = conj(b) * b_dom
-
-        let X = gram_eigh_full(a, (i, k), (j, l))
-            X_cod = replacedimnames(X, "j" => "i", "l" => "k")
-            @test (j, l) ⊆ inds(X)
-            @test X_cod * conj(X) ≈ a
-        end
-
-        let (X, Y) = gram_eigh_full_with_pinv(a, (i, k), (j, l))
-            rank_name = only(setdiff(dimnames(X), ("j", "l")))
-            @test rank_name == only(setdiff(dimnames(Y), ("j", "l")))
-            X_cod = replacedimnames(X, "j" => "i", "l" => "k")
-            @test X_cod * conj(X) ≈ a
-            # Rename one rank dimension so `Y * X` contracts only on
-            # the shared domain names `(j, l)` and leaves a
-            # (rank × rank) named identity.
-            fresh_rank = uniquename(rank_name)
-            X_fresh = replacedimnames(X, rank_name => fresh_rank)
-            YXmat = unname(Y * X_fresh, (rank_name, fresh_rank))
-            @test YXmat ≈ LinearAlgebra.I(size(YXmat, 1))
-        end
-    end
     @testset "tr" begin
         i, j = Index.((2, 3))
         ip, jp = prime(i), prime(j)