diff --git a/Project.toml b/Project.toml index dfafa378..0d7532fc 100644 --- a/Project.toml +++ b/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 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" diff --git a/docs/Project.toml b/docs/Project.toml index fbbbfee7..570b05ff 100644 --- a/docs/Project.toml +++ b/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" diff --git a/src/namedtensoroperator.jl b/src/namedtensoroperator.jl index 92312aca..87591b9a 100644 --- a/src/namedtensoroperator.jl +++ b/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...) diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index 4911616e..38255ed4 100644 --- a/src/tensoralgebra.jl +++ b/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 diff --git a/test/Project.toml b/test/Project.toml index 7e3740f5..912ae1bf 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/test_operator.jl b/test/test_operator.jl index cb377d5d..f1720472 100644 --- a/test/test_operator.jl +++ b/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) diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl index 7c35719c..54b4ddfd 100644 --- a/test/test_tensoralgebra.jl +++ b/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)