From cf431eb322c122156e3587be73e11832d21a6dc8 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Sun, 5 Jul 2026 15:47:38 -0400 Subject: [PATCH 1/8] Add named methods for the Hermitian-safe matrix functions and project_hermitian Adds `AbstractNamedTensor` methods for `sqrth_safe`, `invsqrth_safe`, `sqrth_invsqrth_safe`, and `MatrixAlgebraKit.project_hermitian`, dispatching through TensorAlgebra's labeled forms and returning results with the input's dimension names. Updates the `TensorAlgebra` compat to `0.17`. Co-authored-by: Claude Fable 5 --- Project.toml | 7 +- docs/Project.toml | 2 +- src/abstractnamedtensor.jl | 22 ++-- src/tensoralgebra.jl | 215 ++++++++++++++++++++++++++----------- test/Project.toml | 2 +- test/test_tensoralgebra.jl | 6 +- test/test_tensorkitext.jl | 11 +- 7 files changed, 185 insertions(+), 80 deletions(-) diff --git a/Project.toml b/Project.toml index 0412584..6c4cf0d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ITensorBase" uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" -version = "0.10.8" +version = "0.10.9" authors = ["ITensor developers and contributors"] [workspace] @@ -32,6 +32,9 @@ OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" +[sources.TensorAlgebra] +path = "/Users/mfishman/.julia/dev/TensorAlgebra/.worktrees/mf/project-trailing-axes" + [extensions] ITensorBaseAdaptExt = "Adapt" ITensorBaseMooncakeExt = "Mooncake" @@ -54,7 +57,7 @@ OMEinsumContractionOrders = "1" OrderedCollections = "1.6" Random = "1.10" SimpleTraits = "0.9.4" -TensorAlgebra = "0.16.5" +TensorAlgebra = "0.17" TensorKit = "0.17" TensorOperations = "5.3.1" TermInterface = "2" diff --git a/docs/Project.toml b/docs/Project.toml index 4be8f95..f52464c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,5 +16,5 @@ ITensorBase = "0.10" ITensorFormatter = "0.2.27" Literate = "2" MatrixAlgebraKit = "0.2, 0.3, 0.4, 0.5, 0.6" -TensorAlgebra = "0.16" +TensorAlgebra = "0.17" Test = "1.10" diff --git a/src/abstractnamedtensor.jl b/src/abstractnamedtensor.jl index 7043938..66b8944 100644 --- a/src/abstractnamedtensor.jl +++ b/src/abstractnamedtensor.jl @@ -221,8 +221,10 @@ Base.CartesianIndices(a::AbstractNamedTensor) = CartesianIndices(axes(a)) # Forward `conj` to the underlying so that graded axes flip their sector arrows. # The default `AbstractArray` fallback would broadcast `conj` over elements without # touching the axes, which silently changes the contraction convention for tensors -# with graded (dual-tagged) axes. -Base.conj(a::AbstractNamedTensor) = nameddimsof(a, conj(unnamed(a))) +# with graded (dual-tagged) axes. Routed through `TensorAlgebra.conjugate` (not +# `Base.conj`) so non-`AbstractArray` backends whose native conjugation is +# `adjoint`-shaped (e.g. a `TensorMap`) can overload it without piracy. +Base.conj(a::AbstractNamedTensor) = nameddimsof(a, TensorAlgebra.conjugate(unnamed(a))) # `LinearAlgebra.normalize` infers result eltype via `typeof(first(a)/nrm)`, which # scalar-indexes block-structured storage. `a / norm(a, p)` already preserves names. @@ -399,14 +401,17 @@ end Base.axes(a::AbstractNamedTensor, dimname::Name) = axes(a, dim(a, dimname)) Base.size(a::AbstractNamedTensor, dimname::Name) = size(a, dim(a, dimname)) +# Lowered through `TensorAlgebra.similar_map` (all-codomain, so identical to +# `similar(parent, elt, axes)` for dense) so non-`AbstractArray` backends whose `similar` +# wants a map-shaped space (e.g. a `TensorMap`) allocate through their own overload. function similar_nameddims(a::AbstractNamedTensor, elt::Type, ax) return nameddims( - similar(unnamed(a), elt, unnamed.(Tuple(ax))), + TensorAlgebra.similar_map(unnamed(a), elt, unnamed.(Tuple(ax)), ()), name.(ax) ) end function similar_nameddims(a::AbstractArray, elt::Type, ax) - return nameddims(similar(a, elt, unnamed.(Tuple(ax))), name.(ax)) + return nameddims(TensorAlgebra.similar_map(a, elt, unnamed.(Tuple(ax)), ()), name.(ax)) end # Base.similar gets the eltype at compile time. @@ -715,8 +720,10 @@ function Base.getindex( ) return getindex(a, to_indices(a, (I1, Irest...))...) end +# Routed through `TensorAlgebra.scalar` (default `a[]`) so non-`AbstractArray` backends +# without trivial-sector scalar indexing (e.g. a `TensorMap`) can overload it. function Base.getindex(a::AbstractNamedTensor) - return getindex(unnamed(a)) + return TensorAlgebra.scalar(unnamed(a)) end # Linear indexing. function Base.getindex(a::AbstractNamedTensor, I::Int) @@ -1191,9 +1198,10 @@ end # `sum` is routed to the underlying data rather than left to fall back on the # `mapreduce` method above because some array types (such as graded arrays) define # `Base.sum` directly but not the general `mapreduce`, so the unnamed `sum` is the -# path that works for them. +# path that works for them. Routed through `TensorAlgebra.sum` so non-iterable +# backends (e.g. a `TensorMap`) can overload it. function Base.sum(a::AbstractNamedTensor; kwargs...) - return sum(unnamed(a); kwargs...) + return TensorAlgebra.sum(unnamed(a); kwargs...) end function LinearAlgebra.promote_leaf_eltypes(a::AbstractNamedTensor) diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index 0dfebf3..17b0ae5 100644 --- a/src/tensoralgebra.jl +++ b/src/tensoralgebra.jl @@ -438,6 +438,112 @@ function gram_eigh_full_with_pinv_nameddims( return nameddims(x_unnamed, dimnames_x), nameddims(y_unnamed, dimnames_y) end +""" + TensorAlgebra.MatrixAlgebra.sqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p + +Square root of a named array `a`, interpreting it as an approximately +Hermitian positive semi-definite linear map from the domain to the codomain +dimension names. The result carries the same dimension names as `a`. +Eigenvalues below tolerance are clamped to zero. + +`kwargs` are forwarded to `TensorAlgebra.sqrth_safe` on the underlying +unnamed array (e.g. `atol`, `rtol`). + +See also [`TensorAlgebra.MatrixAlgebra.invsqrth_safe`](@ref) and +[`TensorAlgebra.MatrixAlgebra.sqrth_invsqrth_safe`](@ref). +""" +function MA.sqrth_safe( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + return sqrth_safe_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) +end +function sqrth_safe_nameddims( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + codomain = name.(dimnames_codomain) + domain = name.(dimnames_domain) + p_unnamed = TA.sqrth_safe(unnamed(a), dimnames(a), codomain, domain; kwargs...) + return nameddims(p_unnamed, (codomain..., domain...)) +end + +""" + TensorAlgebra.MatrixAlgebra.invsqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p + +Pseudo-inverse square root of a named array `a`, interpreting it as an +approximately Hermitian positive semi-definite linear map from the domain +to the codomain dimension names. The result carries the same dimension +names as `a`. Eigenvalues below tolerance are clamped to zero +(Moore-Penrose convention). + +`kwargs` are forwarded to `TensorAlgebra.invsqrth_safe` on the underlying +unnamed array (e.g. `atol`, `rtol`). + +See also [`TensorAlgebra.MatrixAlgebra.sqrth_safe`](@ref) and +[`TensorAlgebra.MatrixAlgebra.sqrth_invsqrth_safe`](@ref). +""" +function MA.invsqrth_safe( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + return invsqrth_safe_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) +end +function invsqrth_safe_nameddims( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + codomain = name.(dimnames_codomain) + domain = name.(dimnames_domain) + p_unnamed = TA.invsqrth_safe(unnamed(a), dimnames(a), codomain, domain; kwargs...) + return nameddims(p_unnamed, (codomain..., domain...)) +end + +""" + TensorAlgebra.MatrixAlgebra.sqrth_invsqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p, pinv + +Square root and pseudo-inverse square root of a named array `a` (see +`TensorAlgebra.MatrixAlgebra.sqrth_safe` and +`TensorAlgebra.MatrixAlgebra.invsqrth_safe`), from a single +eigendecomposition. Both results carry the same dimension names as `a`. + +`kwargs` are forwarded to `TensorAlgebra.sqrth_invsqrth_safe` on the underlying +unnamed array (e.g. `atol`, `rtol`). +""" +function MA.sqrth_invsqrth_safe( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + return sqrth_invsqrth_safe_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) +end +function sqrth_invsqrth_safe_nameddims( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + codomain = name.(dimnames_codomain) + domain = name.(dimnames_domain) + p_unnamed, pinv_unnamed = TA.sqrth_invsqrth_safe( + unnamed(a), dimnames(a), codomain, domain; kwargs... + ) + dimnames_p = (codomain..., domain...) + return nameddims(p_unnamed, dimnames_p), nameddims(pinv_unnamed, dimnames_p) +end + +""" + MatrixAlgebraKit.project_hermitian(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> h + +Hermitian part `(m + m') / 2` of a named array `a`, interpreting it as a +linear map `m` from the domain to the codomain dimension names. The result +carries the same dimension names as `a`. +""" +function MAK.project_hermitian( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + return project_hermitian_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) +end +function project_hermitian_nameddims( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + codomain = name.(dimnames_codomain) + domain = name.(dimnames_domain) + h_unnamed = TA.project_hermitian(unnamed(a), dimnames(a), codomain, domain; kwargs...) + return nameddims(h_unnamed, (codomain..., domain...)) +end + """ Base.one(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain) -> Id @@ -513,69 +619,56 @@ end # Projection into a symmetry-restricted named tensor. # -# Shared implementation: strip to `a`'s axes, lower to the TensorAlgebra hook, and reattach the -# names. The dispatch entries below feed this from either a codomain/domain split or a flat list -# (an empty domain). Two split entries per function read the index type from whichever side is -# non-empty (an empty codomain is the all-domain case, the mirror of the empty-domain state), so -# neither an empty codomain nor an empty domain falls through to the unnamed-axis generic. -function project_nameddims(a, codomain_inds, domain_inds) - raw = TA.project(a, unnamed.(codomain_inds), unnamed.(domain_inds)) - return nameddims(raw, (name.(codomain_inds)..., name.(domain_inds)...)) -end -function checked_project_nameddims(a, codomain_inds, domain_inds; kwargs...) - raw = TA.checked_project(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...) - return nameddims(raw, (name.(codomain_inds)..., name.(domain_inds)...)) -end - """ - TensorAlgebra.project(a::AbstractArray, codomain_inds, domain_inds) -> t - TensorAlgebra.project(a::AbstractArray, inds) -> t + TensorAlgebra.project(a::AbstractArray, codomain_inds, domain_inds; kwargs...) -> t + TensorAlgebra.project(a::AbstractArray, inds; kwargs...) -> t Build a named tensor from the dense array `a` by projecting it into the -symmetry-restricted space described by the indices. The three-argument form -takes an explicit codomain/domain split (an operator); the two-argument form -takes a flat list of indices (a state, i.e. an empty domain). The index axes -select the backend: dense ranges give an `Array`, graded ranges a block-sparse -array, and TensorKit spaces a `TensorMap`. `a` is indexed positionally in the -order `(codomain_inds..., domain_inds...)`. - -See `TensorAlgebra.checked_project` for a version that verifies nothing outside -the symmetry-allowed blocks was discarded. +symmetry-restricted space described by the indices, verifying that only a +negligible component of `a` is discarded and throwing an `InexactError` +otherwise (keyword arguments are forwarded to the `isapprox` tolerance +check). The three-argument form takes an explicit codomain/domain split (an +operator); the two-argument form takes a flat list of indices (a state, i.e. +an empty domain). The index axes select the backend: dense ranges give an +`Array`, graded ranges a block-sparse array, and TensorKit spaces a +`TensorMap`. `a` is indexed positionally in the order +`(codomain_inds..., domain_inds...)`. + +`TensorAlgebra.tryproject` returns `nothing` instead of throwing, and +`TensorAlgebra.unchecked_project` skips the verification. """ -function TA.project( - a::AbstractArray, - codomain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}, - domain_inds::Tuple{Vararg{NamedUnitRange}} - ) - return project_nameddims(a, codomain_inds, domain_inds) -end -function TA.project( - a::AbstractArray, - codomain_inds::Tuple{}, - domain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}} - ) - return project_nameddims(a, codomain_inds, domain_inds) -end -function TA.project(a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}) - return project_nameddims(a, inds, ()) -end - -function TA.checked_project( - a::AbstractArray, - codomain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}, - domain_inds::Tuple{Vararg{NamedUnitRange}}; kwargs... - ) - return checked_project_nameddims(a, codomain_inds, domain_inds; kwargs...) -end -function TA.checked_project( - a::AbstractArray, - codomain_inds::Tuple{}, - domain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; kwargs... - ) - return checked_project_nameddims(a, codomain_inds, domain_inds; kwargs...) -end -function TA.checked_project( - a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; kwargs... - ) - return checked_project_nameddims(a, inds, (); kwargs...) +TA.project + +# Strip to `a`'s axes, lower to the TensorAlgebra verb, and reattach the names (passing +# through a `nothing` from `tryproject`). Two split entries per verb read the index type from +# whichever side is non-empty (an empty codomain is the all-domain case, the mirror of the +# empty-domain state), so neither an empty codomain nor an empty domain falls through to the +# unnamed-axis generic; the flat (state) form forwards to the split form. +for f in (:project, :tryproject, :unchecked_project) + @eval begin + function TA.$f( + a::AbstractArray, + codomain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}, + domain_inds::Tuple{Vararg{NamedUnitRange}}; kwargs... + ) + raw = TA.$f(a, unnamed.(codomain_inds), unnamed.(domain_inds); kwargs...) + isnothing(raw) && return nothing + return nameddims(raw, (name.(codomain_inds)..., name.(domain_inds)...)) + end + function TA.$f( + a::AbstractArray, + codomain_inds::Tuple{}, + domain_inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; kwargs... + ) + raw = TA.$f(a, (), unnamed.(domain_inds); kwargs...) + isnothing(raw) && return nothing + return nameddims(raw, name.(domain_inds)) + end + function TA.$f( + a::AbstractArray, inds::Tuple{NamedUnitRange, Vararg{NamedUnitRange}}; + kwargs... + ) + return TA.$f(a, inds, (); kwargs...) + end + end end diff --git a/test/Project.toml b/test/Project.toml index 5a2e453..f8ac67a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -42,7 +42,7 @@ Random = "1.10" SafeTestsets = "0.1" StableRNGs = "1" Suppressor = "0.2" -TensorAlgebra = "0.16" +TensorAlgebra = "0.17" TensorKit = "0.17" TensorOperations = "5.3.1" TermInterface = "2" diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl index 3998e50..5a8f258 100644 --- a/test/test_tensoralgebra.jl +++ b/test/test_tensoralgebra.jl @@ -6,7 +6,7 @@ using MatrixAlgebraKit: left_null, left_orth, left_polar, lq_compact, lq_full, q using StableRNGs: StableRNG using TensorAlgebra.MatrixAlgebra: gram_eigh_full, gram_eigh_full_with_pinv using TensorAlgebra: - TensorAlgebra, checked_project, contract, matricize, project, unmatricize + TensorAlgebra, contract, matricize, project, unchecked_project, unmatricize using Test: @test, @test_broken, @testset @testset "TensorAlgebra (eltype=$(elt))" for elt in @@ -165,8 +165,8 @@ using Test: @test, @test_broken, @testset @test eltype(top) === elt @test Set(dimnames(top)) == Set(name.((prime(i), i))) @test unname(top, (prime(i), i)) == Sz - # `checked_project` accepts the (for dense, always exact) projection - @test unname(checked_project(Sz, (prime(i),), (i,)), (prime(i), i)) == Sz + # `unchecked_project` skips the (for dense, always exact) verification + @test unname(unchecked_project(Sz, (prime(i),), (i,)), (prime(i), i)) == Sz # the two-argument form builds a state (empty domain) v = elt[1, 0] s = project(v, (i,)) diff --git a/test/test_tensorkitext.jl b/test/test_tensorkitext.jl index 24efcd0..6a00c5e 100644 --- a/test/test_tensorkitext.jl +++ b/test/test_tensorkitext.jl @@ -2,7 +2,7 @@ using ITensorBase: ITensorBase, Index, aligndims, dimnames, name, prime, unnamed using LinearAlgebra: norm using MatrixAlgebraKit: qr_compact, svd_compact using StableRNGs: StableRNG -using TensorAlgebra: TensorAlgebra, checked_project, project +using TensorAlgebra: TensorAlgebra, project, unchecked_project using TensorKit: TensorKit, @tensor, AbstractTensorMap, SU2Irrep, U1Irrep, Vect, dim, dual, scalar, space, ←, ⊗ using Test: @test, @test_throws, @testset @@ -127,16 +127,17 @@ using Test: @test, @test_throws, @testset @test space(unnamed(top)) == (W ← W) @test Set(dimnames(top)) == Set(name.((prime(w), w))) - # a charge-breaking operator is projected to zero; `checked_project` rejects the discard - @test norm(unnamed(project(Sx, (prime(w),), (w,)))) == 0 - @test_throws InexactError checked_project(Sx, (prime(w),), (w,); atol = 0, rtol = 0) + # a charge-breaking operator is projected to zero by `unchecked_project`; the checked + # `project` rejects the discard + @test norm(unnamed(unchecked_project(Sx, (prime(w),), (w,)))) == 0 + @test_throws InexactError project(Sx, (prime(w),), (w,); atol = 0, rtol = 0) # the two-argument form builds an all-codomain state; only the trivial-charge component # of the dense vector survives pv = project(elt[1, 0], (w,)) @test unnamed(pv) isa AbstractTensorMap @test norm(unnamed(pv)) ≈ 1 - @test norm(unnamed(project(elt[0, 1], (w,)))) == 0 + @test norm(unnamed(unchecked_project(elt[0, 1], (w,)))) == 0 # the empty-codomain form builds an all-domain `TensorMap` (the mirror case) cobra = project(elt[1, 0], (), (w,)) From 438f522fd94a041cb787f11dba1c56a7093efc81 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Sun, 5 Jul 2026 19:32:58 -0400 Subject: [PATCH 2/8] Generalize the broadcast alignment barrier to non-AbstractArray backends Renames `_broadcast_permuteddims_to` to `_broadcast_permuteddims` and merges it into the existing barrier as one untyped method, taking the rank from `TensorAlgebra.ndims` so operands that are not `AbstractArray`s (such as a `TensorMap`) go through the same path. --- src/broadcast.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/broadcast.jl b/src/broadcast.jl index 14b8b0e..1b65513 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -27,7 +27,7 @@ function broadcasted_unnamed(a::AbstractNamedTensor, names) # common case for the rest) needs no permutation, avoiding a `getperm` allocation and the # identity `permuteddims` wrapper. Skipping it makes a small add several times slower. dimnames(a) == names && return unnamed(a) - return _broadcast_permuteddims_to(unnamed(a), getperm(dimnames(a), names)) + return _broadcast_permuteddims(unnamed(a), getperm(dimnames(a), names)) end # Broadcasting-only alignment: unlike the public `unnamed(a, names)` (which returns a # `Base.PermutedDimsArray`, a full array), this wraps in `TensorAlgebra.PermutedDims`, which stores @@ -35,10 +35,11 @@ end # from the runtime permutation and is a broadcast leaf the linear-combination fold absorbs via # `bipermutedimsopadd!`. `PermutedDims` has almost no array interface, so it stays confined to this # hot path and is never handed back to users. Function barrier: `unnamed(a)` is abstractly typed, -# so dispatching on the concrete array makes `ndims` a compile-time constant for the inferrable -# `ntuple(…, Val(ndims))` permutation. -@noinline function _broadcast_permuteddims_to(array::AbstractArray, perm) - return TA.PermutedDims(array, ntuple(i -> perm[i], Val(ndims(array)))) +# so dispatching on the concrete array makes the rank a compile-time constant for the inferrable +# `ntuple(…, Val(ndims))` permutation. The rank comes from `TensorAlgebra.ndims`, which also +# covers non-`AbstractArray` backends like a `TensorMap`. +@noinline function _broadcast_permuteddims(array, perm) + return TA.PermutedDims(array, ntuple(i -> perm[i], Val(TA.ndims(array)))) end function broadcasted_unnamed(bc::Broadcasted, names) return broadcasted(bc.f, Base.Fix2(broadcasted_unnamed, names).(bc.args)...) From ecb91d2d7260a4a0b67c6920f70361d26748f7b7 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Sun, 5 Jul 2026 21:31:09 -0400 Subject: [PATCH 3/8] Make uniquename mint bare IndexNames uniquename on an IndexName now produces a fresh bare name (new id, no tags, prime level zero) instead of keeping the seed's decoration. A fresh name, such as a factorization's new bond, has no relationship to the seed's tags or prime level, and inheriting them made factorization bond names depend on the input tensor's storage order (the first stored dimension seeded the name). Also documents that the named safe spectral functions require Hermitian input, matching the TensorAlgebra change. Co-authored-by: Claude --- src/index.jl | 7 ++++--- src/tensoralgebra.jl | 22 +++++++++++++--------- src/uniquename.jl | 7 ++++--- 3 files changed, 21 insertions(+), 15 deletions(-) diff --git a/src/index.jl b/src/index.jl index a709c64..73a837d 100644 --- a/src/index.jl +++ b/src/index.jl @@ -33,9 +33,10 @@ function IndexName( ) return IndexName(id, to_tags(tags), plev) end -function uniquename(rng::AbstractRNG, n::IndexName) - return setid(n, uuid4(rng)) -end +# Minting from an existing name deliberately drops its tags and prime level (the generic +# name-instance fallback routes through this type form): a fresh name (a factorization's +# new bond, an operator's fresh codomain) has no relationship to the seed's decoration, +# only to its name type. function uniquename(rng::AbstractRNG, ::Type{<:IndexName}) return IndexName(rng) end diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index 17b0ae5..e0bf230 100644 --- a/src/tensoralgebra.jl +++ b/src/tensoralgebra.jl @@ -441,10 +441,12 @@ end """ TensorAlgebra.MatrixAlgebra.sqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p -Square root of a named array `a`, interpreting it as an approximately -Hermitian positive semi-definite linear map from the domain to the codomain -dimension names. The result carries the same dimension names as `a`. -Eigenvalues below tolerance are clamped to zero. +Square root of a named array `a`, interpreting it as a Hermitian positive +semi-definite linear map from the domain to the codomain dimension names. +The result carries the same dimension names as `a`. Eigenvalues below +tolerance are clamped to zero. The input must be Hermitian: project with +`MatrixAlgebraKit.project_hermitian` first if it is Hermitian only up to +numerical noise. `kwargs` are forwarded to `TensorAlgebra.sqrth_safe` on the underlying unnamed array (e.g. `atol`, `rtol`). @@ -469,11 +471,13 @@ end """ TensorAlgebra.MatrixAlgebra.invsqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p -Pseudo-inverse square root of a named array `a`, interpreting it as an -approximately Hermitian positive semi-definite linear map from the domain -to the codomain dimension names. The result carries the same dimension -names as `a`. Eigenvalues below tolerance are clamped to zero -(Moore-Penrose convention). +Pseudo-inverse square root of a named array `a`, interpreting it as a +Hermitian positive semi-definite linear map from the domain to the codomain +dimension names. The result carries the same dimension names as `a`. +Eigenvalues below tolerance are clamped to zero (Moore-Penrose convention). +The input must be Hermitian: project with +`MatrixAlgebraKit.project_hermitian` first if it is Hermitian only up to +numerical noise. `kwargs` are forwarded to `TensorAlgebra.invsqrth_safe` on the underlying unnamed array (e.g. `atol`, `rtol`). diff --git a/src/uniquename.jl b/src/uniquename.jl index e536196..7da0458 100644 --- a/src/uniquename.jl +++ b/src/uniquename.jl @@ -6,9 +6,10 @@ using Random: AbstractRNG, RandomDevice Mint a fresh, unique name. Given an existing `name` (or a name `type`), produce a new name of the same flavor that is distinct from any other, for example to label a freshly -generated dimension in a matrix factorization. Randomness defaults to OS entropy -(`Random.RandomDevice`) so that minting a name neither perturbs nor is perturbed by the -numerical RNG. Pass an explicit `rng` for a reproducible name. +generated dimension in a matrix factorization. For decorated name types the fresh name +is bare: an `IndexName` seed's tags and prime level are not inherited. Randomness +defaults to OS entropy (`Random.RandomDevice`) so that minting a name neither perturbs +nor is perturbed by the numerical RNG. Pass an explicit `rng` for a reproducible name. # Examples From 8b8d9285f914583225114f6ff8095ac3f07f38ef Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 6 Jul 2026 09:45:30 -0400 Subject: [PATCH 4/8] Collapse the named safe-spectral wrappers into one loop The named `sqrth_safe`, `invsqrth_safe`, and `project_hermitian` methods for `AbstractNamedTensor` share the same body (lower to the unnamed array, reattach names codomain-first), so generate them from one `@eval` loop. `sqrth_invsqrth_safe` stays a standalone method since it fans the names over a returned pair. The docstrings remain standalone binding docs. Co-authored-by: Claude --- src/tensoralgebra.jl | 72 ++++++++++++++------------------------------ 1 file changed, 23 insertions(+), 49 deletions(-) diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index e0bf230..de83c24 100644 --- a/src/tensoralgebra.jl +++ b/src/tensoralgebra.jl @@ -454,19 +454,7 @@ unnamed array (e.g. `atol`, `rtol`). See also [`TensorAlgebra.MatrixAlgebra.invsqrth_safe`](@ref) and [`TensorAlgebra.MatrixAlgebra.sqrth_invsqrth_safe`](@ref). """ -function MA.sqrth_safe( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - return sqrth_safe_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) -end -function sqrth_safe_nameddims( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - codomain = name.(dimnames_codomain) - domain = name.(dimnames_domain) - p_unnamed = TA.sqrth_safe(unnamed(a), dimnames(a), codomain, domain; kwargs...) - return nameddims(p_unnamed, (codomain..., domain...)) -end +MA.sqrth_safe """ TensorAlgebra.MatrixAlgebra.invsqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p @@ -485,19 +473,7 @@ unnamed array (e.g. `atol`, `rtol`). See also [`TensorAlgebra.MatrixAlgebra.sqrth_safe`](@ref) and [`TensorAlgebra.MatrixAlgebra.sqrth_invsqrth_safe`](@ref). """ -function MA.invsqrth_safe( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - return invsqrth_safe_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) -end -function invsqrth_safe_nameddims( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - codomain = name.(dimnames_codomain) - domain = name.(dimnames_domain) - p_unnamed = TA.invsqrth_safe(unnamed(a), dimnames(a), codomain, domain; kwargs...) - return nameddims(p_unnamed, (codomain..., domain...)) -end +MA.invsqrth_safe """ TensorAlgebra.MatrixAlgebra.sqrth_invsqrth_safe(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> p, pinv @@ -510,22 +486,7 @@ eigendecomposition. Both results carry the same dimension names as `a`. `kwargs` are forwarded to `TensorAlgebra.sqrth_invsqrth_safe` on the underlying unnamed array (e.g. `atol`, `rtol`). """ -function MA.sqrth_invsqrth_safe( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - return sqrth_invsqrth_safe_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) -end -function sqrth_invsqrth_safe_nameddims( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - codomain = name.(dimnames_codomain) - domain = name.(dimnames_domain) - p_unnamed, pinv_unnamed = TA.sqrth_invsqrth_safe( - unnamed(a), dimnames(a), codomain, domain; kwargs... - ) - dimnames_p = (codomain..., domain...) - return nameddims(p_unnamed, dimnames_p), nameddims(pinv_unnamed, dimnames_p) -end +MA.sqrth_invsqrth_safe """ MatrixAlgebraKit.project_hermitian(a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs...) -> h @@ -534,18 +495,31 @@ Hermitian part `(m + m') / 2` of a named array `a`, interpreting it as a linear map `m` from the domain to the codomain dimension names. The result carries the same dimension names as `a`. """ -function MAK.project_hermitian( - a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... - ) - return project_hermitian_nameddims(a, dimnames_codomain, dimnames_domain; kwargs...) +MAK.project_hermitian + +# The named forms above: lower to the corresponding tensor-level TensorAlgebra function on +# the unnamed array and reattach the names, codomain first. `sqrth_invsqrth_safe` differs +# only in fanning the names out over its result pair. +for (M, f) in ((MA, :sqrth_safe), (MA, :invsqrth_safe), (MAK, :project_hermitian)) + @eval function $M.$f( + a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... + ) + codomain = name.(dimnames_codomain) + domain = name.(dimnames_domain) + p_unnamed = TA.$f(unnamed(a), dimnames(a), codomain, domain; kwargs...) + return nameddims(p_unnamed, (codomain..., domain...)) + end end -function project_hermitian_nameddims( +function MA.sqrth_invsqrth_safe( a::AbstractNamedTensor, dimnames_codomain, dimnames_domain; kwargs... ) codomain = name.(dimnames_codomain) domain = name.(dimnames_domain) - h_unnamed = TA.project_hermitian(unnamed(a), dimnames(a), codomain, domain; kwargs...) - return nameddims(h_unnamed, (codomain..., domain...)) + p_unnamed, pinv_unnamed = TA.sqrth_invsqrth_safe( + unnamed(a), dimnames(a), codomain, domain; kwargs... + ) + dimnames_p = (codomain..., domain...) + return nameddims(p_unnamed, dimnames_p), nameddims(pinv_unnamed, dimnames_p) end """ From d65fd221a73425d0da9d85a248b9382a391d9901 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 6 Jul 2026 10:53:05 -0400 Subject: [PATCH 5/8] Pin the TensorAlgebra branch while the 0.17 support is in review Co-authored-by: Claude Fable 5 --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6c4cf0d..be4ee1c 100644 --- a/Project.toml +++ b/Project.toml @@ -33,7 +33,8 @@ TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" [sources.TensorAlgebra] -path = "/Users/mfishman/.julia/dev/TensorAlgebra/.worktrees/mf/project-trailing-axes" +rev = "mf/project-trailing-axes" +url = "https://github.com/ITensor/TensorAlgebra.jl" [extensions] ITensorBaseAdaptExt = "Adapt" From 376b8e8d8536f19d073eafe060882b4e5611edd1 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 6 Jul 2026 13:12:24 -0400 Subject: [PATCH 6/8] Keep tags and prime level in uniquename on an existing name `uniquename` on an existing `IndexName` now keeps its tags and prime level, minting only a fresh id (the legacy `sim`). Pass the name type instead of an instance to mint a bare name with no tags and prime level zero. The factorization bond names pass the type, since a fresh bond has no relationship to the input's decoration. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/index.jl | 11 +++++++---- src/tensoralgebra.jl | 18 +++++++++--------- src/uniquename.jl | 14 ++++++++------ test/test_basics.jl | 23 +++++++++++++++++++++-- 4 files changed, 45 insertions(+), 21 deletions(-) diff --git a/src/index.jl b/src/index.jl index 73a837d..4f31493 100644 --- a/src/index.jl +++ b/src/index.jl @@ -33,10 +33,13 @@ function IndexName( ) return IndexName(id, to_tags(tags), plev) end -# Minting from an existing name deliberately drops its tags and prime level (the generic -# name-instance fallback routes through this type form): a fresh name (a factorization's -# new bond, an operator's fresh codomain) has no relationship to the seed's decoration, -# only to its name type. +# `uniquename` on an existing `IndexName` keeps its tags and prime level, minting only a +# fresh id (the legacy `sim`). The type form drops them: a factorization bond or a fresh +# operator leg has no relationship to any seed's decoration, so its callers pass the name +# type to opt out of inheriting it. +function uniquename(rng::AbstractRNG, name::IndexName) + return IndexName(rng; tags = tags_stored(name), plev = plev(name)) +end function uniquename(rng::AbstractRNG, ::Type{<:IndexName}) return IndexName(rng) end diff --git a/src/tensoralgebra.jl b/src/tensoralgebra.jl index de83c24..bcc239a 100644 --- a/src/tensoralgebra.jl +++ b/src/tensoralgebra.jl @@ -136,7 +136,7 @@ for f in [ domain = name.(dimnames_domain) x_unnamed, y_unnamed = TA.$f(unnamed(a), dimnames(a), codomain, domain; kwargs...) - name_x = uniquename(dimnames(a, 1)) + name_x = uniquename(dimnametype(a)) name_y = name_x dimnames_x = (codomain..., name_x) dimnames_y = (name_y, domain...) @@ -175,8 +175,8 @@ for f in [:svd_compact, :svd_full, :svd_trunc] u_unnamed, s_unnamed, v_unnamed = TA.$f( unnamed(a), dimnames(a), codomain, domain; kwargs... ) - name_u = uniquename(dimnames(a, 1)) - name_v = uniquename(dimnames(a, 1)) + name_u = uniquename(dimnametype(a)) + name_v = uniquename(dimnametype(a)) dimnames_u = (codomain..., name_u) dimnames_s = (name_u, name_v) dimnames_v = (name_v, domain...) @@ -249,8 +249,8 @@ for f in [:eigh_full, :eig_full, :eigh_trunc, :eig_trunc] d_unnamed, v_unnamed = TA.$f( unnamed(a), dimnames(a), codomain, domain; kwargs... ) - name_d = uniquename(dimnames(a, 1)) - name_d′ = uniquename(name_d) + name_d = uniquename(dimnametype(a)) + name_d′ = uniquename(dimnametype(a)) name_v = name_d dimnames_d = (name_d′, name_d) dimnames_v = (domain..., name_v) @@ -294,7 +294,7 @@ function left_null_nameddims( codomain = name.(dimnames_codomain) domain = name.(dimnames_domain) n_unnamed = TA.left_null(unnamed(a), dimnames(a), codomain, domain; kwargs...) - name_n = uniquename(dimnames(a, 1)) + name_n = uniquename(dimnametype(a)) dimnames_n = (codomain..., name_n) return nameddims(n_unnamed, dimnames_n) end @@ -319,7 +319,7 @@ function right_null_nameddims( codomain = name.(dimnames_codomain) domain = name.(dimnames_domain) n_unnamed = TA.right_null(unnamed(a), dimnames(a), codomain, domain; kwargs...) - name_n = uniquename(dimnames(a, 1)) + name_n = uniquename(dimnametype(a)) dimnames_n = (name_n, domain...) return nameddims(n_unnamed, dimnames_n) end @@ -376,7 +376,7 @@ function gram_eigh_full_nameddims( codomain = name.(dimnames_codomain) domain = name.(dimnames_domain) x_unnamed = TA.gram_eigh_full(unnamed(a), dimnames(a), codomain, domain; kwargs...) - name_x = uniquename(dimnames(a, 1)) + name_x = uniquename(dimnametype(a)) dimnames_x = (domain..., name_x) return nameddims(x_unnamed, dimnames_x) end @@ -432,7 +432,7 @@ function gram_eigh_full_with_pinv_nameddims( x_unnamed, y_unnamed = TA.gram_eigh_full_with_pinv( unnamed(a), dimnames(a), codomain, domain; kwargs... ) - name_xy = uniquename(dimnames(a, 1)) + 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) diff --git a/src/uniquename.jl b/src/uniquename.jl index 7da0458..ba9a555 100644 --- a/src/uniquename.jl +++ b/src/uniquename.jl @@ -4,12 +4,14 @@ using Random: AbstractRNG, RandomDevice uniquename([rng,] name) uniquename([rng,] type::Type) -Mint a fresh, unique name. Given an existing `name` (or a name `type`), produce a new -name of the same flavor that is distinct from any other, for example to label a freshly -generated dimension in a matrix factorization. For decorated name types the fresh name -is bare: an `IndexName` seed's tags and prime level are not inherited. Randomness -defaults to OS entropy (`Random.RandomDevice`) so that minting a name neither perturbs -nor is perturbed by the numerical RNG. Pass an explicit `rng` for a reproducible name. +Mint a fresh, unique name. Given an existing `name`, produce a new name of the same flavor +that is distinct from any other, for example to label a freshly generated dimension in a +matrix factorization. Decoration carried by the seed is kept: `uniquename` on an `IndexName` +keeps its tags and prime level, minting only a fresh id. Pass the name `type` instead of an +instance to mint a bare name (for `IndexName`, no tags and prime level zero), for a fresh +dimension that should not inherit any seed's decoration. Randomness defaults to OS entropy +(`Random.RandomDevice`) so that minting a name neither perturbs nor is perturbed by the +numerical RNG. Pass an explicit `rng` for a reproducible name. # Examples diff --git a/test/test_basics.jl b/test/test_basics.jl index a4122a1..748d3ee 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -1,6 +1,6 @@ using ITensorBase: ITensorBase, Index, IndexName, NamedTensor, dimnametype, gettag, hastag, - id, inds, mapinds, name, named, plev, prime, setplev, settag, tags, unname, unnamed, - unsettag + id, inds, mapinds, name, named, plev, prime, setplev, settag, tags, uniquename, unname, + unnamed, unsettag using Test: @test, @test_broken, @test_throws, @testset using UUIDs: UUID @@ -47,6 +47,25 @@ using UUIDs: UUID @test gettag(n, "X") isa AbstractString @test ITensorBase.tags_stored(n)[:X] === :Y end + @testset "uniquename" begin + i = settag(prime(Index(2)), "X", "Y") + # On an instance, only the id is fresh: tags and prime level are kept. + n = uniquename(name(i)) + @test n != name(i) + @test id(n) != id(name(i)) + @test plev(n) == plev(i) + @test tags(n) == tags(i) + # On the name type, a bare name: no tags, prime level zero. + m = uniquename(IndexName) + @test m isa IndexName + @test plev(m) == 0 + @test isempty(tags(m)) + # On an `Index`, recurse into the name and keep its decoration. + i′ = uniquename(i) + @test name(i′) != name(i) + @test plev(i′) == plev(i) + @test tags(i′) == tags(i) + end @testset "Index basics" begin i = Index(2) @test plev(i) == 0 From 42f94af4fc8f01f7247fb054bf018ce35ccf9242 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Mon, 6 Jul 2026 20:21:12 -0400 Subject: [PATCH 7/8] Conjugate named tensors through broadcasting Define `conj` on an `AbstractNamedTensor` as `conj.(a)`, so eager conjugation reuses the linear-combination broadcast machinery instead of a separate backend hook. `conj.` lowers to a `ConjBroadcasted` leaf (introduced in https://github.com/ITensor/TensorAlgebra.jl/pull/204) that the op primitives absorb, materializing into an all-codomain destination with dualized axes, so a named tensor wrapping a `TensorMap` can be conjugated. The broadcast destination is now allocated from the flattened expression's axes rather than the prototype's, since an axis-dualizing operand makes them differ. --- src/abstractnamedtensor.jl | 15 ++++++++------- src/broadcast.jl | 9 ++++++++- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/abstractnamedtensor.jl b/src/abstractnamedtensor.jl index 66b8944..7e01bc9 100644 --- a/src/abstractnamedtensor.jl +++ b/src/abstractnamedtensor.jl @@ -218,13 +218,14 @@ Base.zero(a::AbstractNamedTensor) = nameddimsof(a, zero(unnamed(a))) # `AbstractArray` fallback did through `axes`). Base.CartesianIndices(a::AbstractNamedTensor) = CartesianIndices(axes(a)) -# Forward `conj` to the underlying so that graded axes flip their sector arrows. -# The default `AbstractArray` fallback would broadcast `conj` over elements without -# touching the axes, which silently changes the contraction convention for tensors -# with graded (dual-tagged) axes. Routed through `TensorAlgebra.conjugate` (not -# `Base.conj`) so non-`AbstractArray` backends whose native conjugation is -# `adjoint`-shaped (e.g. a `TensorMap`) can overload it without piracy. -Base.conj(a::AbstractNamedTensor) = nameddimsof(a, TensorAlgebra.conjugate(unnamed(a))) +# Eager conjugation reuses the linear-combination broadcast machinery: `conj.(a)` lowers to a +# `ConjBroadcasted` leaf that the op primitives absorb (folding `conj` into their `op` flag) and +# materializes into an all-codomain destination with dualized axes. This keeps one conjugation +# mechanism shared by the lazy `conj.` path and eager `conj`, so a non-`AbstractArray` backend +# (a `TensorMap`) needs no separate eager hook. The default `AbstractArray` `conj` would instead +# map over elements without dualizing the axes, silently changing the contraction convention for +# graded tensors. +Base.conj(a::AbstractNamedTensor) = conj.(a) # `LinearAlgebra.normalize` infers result eltype via `typeof(first(a)/nrm)`, which # scalar-indexes block-structured storage. `a / norm(a, p)` already preserves names. diff --git a/src/broadcast.jl b/src/broadcast.jl index 1b65513..e98a19f 100644 --- a/src/broadcast.jl +++ b/src/broadcast.jl @@ -86,10 +86,17 @@ end # call re-specializes on the concrete runtime types and everything below is type-stable # (`eltype(lb)` is now inferrable, no runtime `promote_op`). Inlining the body into `copy` # instead costs one extra allocation per call. +# +# Allocate from `axes(lb)`, the flattened expression's axes, rather than the prototype's own: +# an axis-changing operand (a `conj` leaf dualizes its axes) makes them differ, and the +# destination must match the expression. All axes go in the codomain (empty domain), the +# all-codomain output convention `@tensor` uses for an unbipartitioned left-hand side; on a +# non-bipartitioned backend (a dense array) `similar_map` with an empty domain is a plain +# `similar` over `axes(lb)`. function _copy_unnamed(bc_unnamed, prototype) lb = TA.tryflattenlinear(bc_unnamed) isnothing(lb) && return copy(bc_unnamed) - return copyto!(similar(prototype, eltype(lb)), lb) + return copyto!(TA.similar_map(prototype, eltype(lb), axes(lb), ()), lb) end # `Base.Broadcast.materialize!` otherwise reconstructs the broadcast over `axes(dest)` and From efaca1085558ea73ebe68848c2b010bfebdb548b Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Tue, 7 Jul 2026 15:57:30 -0400 Subject: [PATCH 8/8] Remove [sources] pin on TensorAlgebra now that 0.17 is registered --- Project.toml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Project.toml b/Project.toml index be4ee1c..0d15309 100644 --- a/Project.toml +++ b/Project.toml @@ -32,10 +32,6 @@ OMEinsumContractionOrders = "6f22d1fd-8eed-4bb7-9776-e7d684900715" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" -[sources.TensorAlgebra] -rev = "mf/project-trailing-axes" -url = "https://github.com/ITensor/TensorAlgebra.jl" - [extensions] ITensorBaseAdaptExt = "Adapt" ITensorBaseMooncakeExt = "Mooncake"