Block-aware methods on graded storage#173
Conversation
TensorKitSectors defines the dual representation via `Base.conj` on every sector type (and `dual(::Sector) = conj(a)` is an alias of it), and never defines `Base.adjoint` on sectors. The `Base.adjoint(::SectorRange) = dual` overload in GradedArrays drifted away from that convention, and the `U1(0)'`, `SU2(1//2)'` shorthand for "dual sector" overloaded `'` (adjoint) onto irrep labels, which is not a mathematically meaningful operation on labels. Drops `Base.adjoint = dual` on `GradedOneTo`, `SectorOneTo`, and `SectorRange`. Adds `Base.conj = dual` on the same three. Migrates every `U1(0)'` / `SU2(1//2)'` test call site (and the in-source docstring example) to `conj(U1(0))` / `conj(SU2(1//2))`. Pairs with NamedDimsArrays.jl PR #229's `Base.conj(::AbstractNamedUnitRange)` forwarder, so a named graded axis's `conj` propagates through the wrapper and flips sector arrows on the underlying axis. This is what `TensorAlgebra.similar_map` needs to compose with GradedArrays-backed prototypes.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #173 +/- ##
==========================================
+ Coverage 81.95% 82.90% +0.95%
==========================================
Files 22 22
Lines 1607 1767 +160
==========================================
+ Hits 1317 1465 +148
- Misses 290 302 +12
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…d `+`-mapreduce on `AbelianGradedArray`
The default `Base`/`LinearAlgebra` methods on `AbstractArray` iterate or scalar-index, both of which `AbstractGradedArray` rejects. This patch fills in the block-wise specializations needed for the common ops on `AbelianGradedArray`: `copy`, `conj` (which also flips axis duality, matching the existing `Base.conj` on the axis types), `norm`, `dot`, `:*`/`:/` with a scalar, and `mapreduce` over `+`/`add_sum` (so `sum` works without iterating). Unstored blocks are treated as implicit zero, which is correct for `+`-mapreduce, `norm`, and `dot`.
Also: `Base.fill!(::FusedGradedVector, val)` now fills every stored block with `val`. The previous `iszero(val) || throw(...)` was too strict, since each block of a fused vector is a per-sector segment whose entries are independently meaningful (`val * I` in the sector-fused basis is the canonical use). `Base.fill!(::AbelianGradedArray, val)` keeps the `iszero` guard, because filling stored blocks with a nonzero value produces a tensor that mixes the value on symmetry-allowed positions with implicit zeros on symmetry-forbidden ones, which is not a meaningful object.
And: scalar `TensorAlgebra.unmatricize(::SectorFusion, ::FusedGradedMatrix, (), ())` is implemented (was an explicit `error`). Reads the trivial-sector entry of `m` into a 0-D `Array{T}`, which is what the unfused codomain/domain axes describe.
Co-authored-by: Claude <noreply@anthropic.com>
…tors on `AbelianGradedArray`
`TensorAlgebra.projectto!(dest::AbelianGradedArray, src::AbstractArray)` copies the symmetry-allowed entries of `src` into `dest` block-by-block and zeroes anything else. Magnitude-blind by contract: any weight on forbidden blocks is dropped silently, leaving the tolerance check to `checked_projectto!`.
`Base.Array(a::AbelianGradedArray)` materializes the graded array into a dense `Array` by laying stored blocks into a zero-initialized destination. This is the densification primitive the default `checked_projectto!` reconstruct-compare check needs.
Two `similar(::AbstractArray|::StridedArray, ::Type{T}, ::Tuple{GradedOneTo, ...})` methods make `similar_map(prototype, T, codomain_axes, domain_axes)` route to an `AbelianGradedArray` allocation whenever the axes are graded, regardless of the prototype's storage type. Without these, dispatch is ambiguous against `BlockArrays`' `StridedArray`-specific overloads.
Co-authored-by: Claude <noreply@anthropic.com>
…path - `LinearAlgebra.isdiag` on `AbelianGradedMatrix` and `FusedGradedMatrix` - `Base.mapreduce` on `FusedGradedVector` - `Base.map` on `FusedGradedVector` - `MatrixAlgebraKit.diagonal` on `FusedGradedVector` - Base matrix functions (`sqrt`, `exp`, `log`, ...) on `FusedGradedMatrix` - `TensorAlgebra.checked_projectto!` on `AbelianGradedArray` Co-authored-by: Claude <noreply@anthropic.com>
- Add `const AbelianGradedMatrix = AbelianGradedArray{T,2,D,S}` for the block-aware `LinearAlgebra.isdiag` dispatch that uses it (precompile was failing without the alias defined at parse time).
- Define `LinearAlgebra.dot(::FusedGradedMatrix, ::FusedGradedMatrix)` block-wise as the sum of `length(c) * dot(A.blocks[c], B.blocks[c])` over shared sectors. The `length(c)` factor matches the per-sector multiplicity convention in `LinearAlgebra.norm`, so `dot(A, A) == norm(A, 2)^2`.
For an `AbstractGradedArray{T, N}`, the splatted method `view(::AbstractGradedArray{T, N}, ::Vararg{Block{1}, N})` collapses to a single `Block{1}` argument when N=1, which is ambiguous with the single-`Block{N}` method. Add explicit `Block{1}`-only methods on `AbstractGradedArray{T, 1}` and `AbelianGradedArray{T, 1}` for `view`, `getindex`, and `setindex!` to resolve the dispatch.
Co-authored-by: Claude <noreply@anthropic.com>
Iterate `values(a.blockdata)` directly instead of going through `eachblockstoredindex(a)` + `view`, which routed through the AbstractGradedArray-generic path and was the wrong shape for the new top-level array ops. `fill!` now actually fills with any value block-wise rather than rejecting nonzero arguments. Co-authored-by: Claude <noreply@anthropic.com>
`FusedGradedMatrix` inherits `BroadcastStyle = GradedStyle{2}()` from `AbstractGradedArray`, but the generic `bipermutedimsopadd!(::AbstractGradedArray, ...)` path on which the broadcast machinery depends assumes the cartesian `(row_block, col_block)` block layout of `AbelianGradedArray`. On a `FusedGradedMatrix`, whose blocks live in a `Dictionary{Sector, Matrix}` keyed at coupled sectors, that path writes to the wrong slots and silently produces wrong results. `A - B` on the residual of an SVD reconstruction was the symptom that surfaced this.
Define `+` and `-` directly block-wise on the union of stored sectors, so a sector present on one side and absent on the other is treated as a zero block. The proper fix is to make broadcasting itself work on `FusedGradedMatrix` (either a per-method specialization of `bipermutedimsopadd!` or a dedicated `FusedGradedStyle`), tracked as a follow-up.
Co-authored-by: Claude <noreply@anthropic.com>
`MAK.truncate(svd_trunc!, ::NTuple{3, FusedGradedMatrix}, ...)` previously kept every coupled sector in the bond, including those whose singular values were all truncated to zero. Downstream gate applications that fused against this bond then produced row/col block indices outside the destination's `allowedblocks` and failed deep in `unmatricize` as `Block (i, j, k) is not stored`.
Drop sectors with zero retained columns from the bond axis only. For U the row (codomain) axis is the input codomain and keeps its full sector set, and only the column (bond) axis shrinks. Analogously for Vᴴ. For S, both sides are the bond and shrink together. Matches the `truncate_space` convention used by `TensorMap` factorizations.
Co-authored-by: Claude <noreply@anthropic.com>
Match the project convention of `using Module: Module [as Alias]` over `import Module [as Alias]`. No behavior change. Co-authored-by: Claude <noreply@anthropic.com>
GradedArrays now references `TensorAlgebra.projectto!`, added in ITensor/TensorAlgebra.jl#177. Without the pin, CI resolves against the registered TensorAlgebra 0.9.4 which does not have `projectto!`, and precompilation fails. Drop the pin once that branch registers. Co-authored-by: Claude <noreply@anthropic.com>
The default `copyto!(::AbstractArray, ::AbstractArray)` fallback uses scalar indexing, which isn't supported on `AbstractGradedArray`. With this specialization, an `AbelianGradedArray → AbelianGradedArray` copy goes through the block dictionary directly. Surfaces in broadcasts that collapse to the bare source array, e.g. `conj.(a)` on a `Real`-eltype `AbelianGradedArray`: `tryflattenlinear` returns the source itself (since `conj` is a no-op on reals), and the subsequent `copyto!(dest, src)` then needed this method. Co-authored-by: Claude <noreply@anthropic.com>
Reject all broadcasts on `FusedGradedMatrix` with an explicit error via a new `FusedGradedStyle`. The generic `GradedStyle` path lowers to `bipermutedimsopadd!` over a cartesian block layout, but `FusedGradedMatrix` stores blocks keyed by coupled sector, so the result was silently wrong. The error points at `+`/`-` or block-wise ops as the supported alternatives until a sector-keyed broadcast path lands. Add block-wise `Base.:*(::FusedGradedMatrix, ::Number)`, the reverse, and `Base.:/` so the common scalar-multiply cases keep working without going through broadcast. Co-authored-by: Claude <noreply@anthropic.com>
Overloads `TensorAlgebra.trivialrange(::Type{<:GradedOneTo})` to return a charge-0 one-dimensional graded range, the identity for `tensor_product_axis` on graded ranges, by delegating to the existing `trivial(R)` method.
Co-authored-by: Claude <noreply@anthropic.com>
Dispatch to each block's `rand!` / `randn!`, alongside the existing block-aware `fill!`, `zero!`, and `scale!`. The generic `AbstractArray` fallbacks go through scalar indexing, which errors on graded arrays. Co-authored-by: Claude <noreply@anthropic.com>
Commit 41a6594 changed `Base.fill!(::AbelianGradedArray, val)` to fill stored blocks block-wise with any value rather than throwing on nonzero. The test still expected an `ArgumentError`. Co-authored-by: Claude <noreply@anthropic.com>
`mortar_axis(::AbstractVector{SectorOneTo{S}})` and `mortar_axis(::AbstractVector{GradedOneTo{S}})` now return `GradedOneTo(S[], Int[])` on empty input. The signatures are parameterized on the sector type `S` so it can be recovered from the eltype.
`sqrt(::FusedGradedMatrix)` (and the rest of the `MATRIX_FUNCTIONS` loop) passes `init = eltype(A)` to the eltype-promoting `mapreduce`, so the reduction works when the block list is empty.
`Base.isequal(::GradedOneTo, ::GradedOneTo)` and `Base.hash` treat two empty graded ranges as equal regardless of `isdual`, since with no sectors there is nothing to flip.
These come up in a U(1) gate-application QR whose codomain index list is empty, producing an empty bond axis that flows through `tensor_product_gradedrange`, `unmatricize`, and `check_input(*, FGM, FGM)`.
`contractopadd!(::Matricize{SectorFusion}, ...)` now verifies that every contracted-axis pair satisfies `dual(axes(a1)[i]) == axes(a2)[j]`. A same-`isdual` pairing (or a sector-content mismatch) was previously aligned block-wise without complaint, producing a number that did not correspond to the intended contraction.
Skipped when the sector's `BraidingStyle` is fermionic, since the supertrace formalism intentionally pairs same-`isdual` axes together with `contraction_twist!` to pick up the right phase.
Same behavior as the previous commit's freestanding helper, but routed through `TensorAlgebra.check_input(::typeof(contract), ::AbstractGradedArray, ...)` so it follows the existing GradedArrays convention for graded-specific input validation. `check_input(contract!)` already delegates to `check_input(contract)`, so the explicit call in `contractopadd!` is no longer needed.
… too Factor the per-axis-pair logic into `axes_match_for_contraction(ax1, ax2)` and call it from `check_input(::typeof(contract), ::AbstractGradedArray, ...)`. The helper always requires sector content to match modulo `isdual`. For bosonic braiding it additionally requires opposite `isdual`, while for fermionic braiding `ax1 == ax2` is also accepted to leave room for the supertrace formalism's same-`isdual` pairings.
Adds graded-array round-trip tests for `TensorAlgebra.svd` and `gram_eigh_full_with_pinv`. The matrix-level tests in TensorAlgebra missed the unmatricize-axis bugs the TA branch fixed (`axes_S` taken from U/Vᴴ rather than S itself, and `axes_Y` not conjugated), because the bugs only surface when the codomain and domain axes are conj-related rather than identical, which only the graded path exercises. Both round-trips are routed through `contract`. The natural `U * S * Vᴴ` and `X * X'` chains fall into the generic scalar-indexing matmul on `AbstractGradedMatrix` and are `@test_broken` pending a block-aware `*`. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
…dArray`
Adds `Base.rand`/`Base.randn(rng, T, ax)` constructors for
`Tuple{GradedOneTo, ...}` axes, a 3-arg
`Random.rand!(rng, ::AbelianGradedArray, ::Random.Sampler)`, and `Base.iszero`.
The constructors allocate an `AbelianGradedArray{T}` with the right block
structure and fill via the in-place methods. The 3-arg `rand!` replaces the
prior 2-arg form so that the `rand!(A)` / `rand!(A, X)` / `rand!(rng, A, X)`
shims, which all converge on the 3-arg `Sampler` path, reach the block-aware
fill. The generic `AbstractArray` paths for all four scalar-index and throw.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Add `FI.permuteddims(a::T, perm) = permutedims(a, perm)` on `AbstractGradedArray`, `AbstractSectorArray`, and `AbstractSectorDelta`. The default `FI.permuteddims` returns a `Base.PermutedDimsArray` view, which falls through GradedArrays' block-aware `bipermutedimsopadd!` overloads (the wrapper itself is not a graded array) and reaches TensorAlgebra's generic implementation. That generic implementation broadcasts back into the graded destination, which lowers through the `GradedStyle` pipeline again and infinite-recurses. NamedDimsArrays' broadcast alignment uses `FI.permuteddims` to align same-name-different-order operands, so any such broadcast on an `AbstractGradedArray` would hang without an eager override. Defining the override at the abstract types covers every concrete GradedArrays array type at once and lets the previous narrower per-subtype overrides go away. Handling lazy permutations systematically, and the future of the `FI.permuteddims` interface itself, remain as follow-ups.
With the GradedArrays-side broadcast hang fixed (ITensor/GradedArrays.jl#173), both legs of the symmetry-looped smoke test now use `prod(gated) ≈ NDA.apply(gate, prod(state))` directly. The `exact_apply_check` helper that special-cased the `:u1` leg with a weaker shape / vertex-set invariant is removed.
Cover behaviors that previously had no direct tests: the `contract` input check that rejects contracted-axis pairs with mismatched duality, `projectto!` / `checked_projectto!` on `AbelianGradedArray`, the `conj` axis-duality flip, `isdiag` on `AbelianGradedMatrix`, and dropping fully truncated sectors from the bond in truncated SVD.
The `Block{N}` view method and its `Block{1}` disambiguator duplicated the `haskey`/sectors/wrap logic, so both now delegate to a shared `view_abelian(a, bk)`.
Define `scale!`, `zero!`, and `fill!` once on `AbstractGradedArray` via the `eachblockstoredindex`/`view` block interface, replacing the per-type versions on `AbelianGradedArray`, `FusedGradedMatrix`, and `FusedGradedVector`. `fill!` is now uniformly permissive: it fills the stored, symmetry-allowed blocks with any value. `AbelianGradedArray` scalar `*` and `/` route through the Base `AbstractArray`-scalar broadcasting, so the dedicated overrides are removed. `FusedGradedMatrix` keeps explicit scalar methods because its broadcasting is disabled, with a TODO to drop them once that is supported. `copyto!` between `AbelianGradedArray`s with matching axes copies each block into the existing destination buffer instead of clearing and reallocating, since matching axes imply matching allocated blocks. `checked_projectto!` forwards its keyword arguments to `isapprox`.
Check axes equality in `dot` on `AbelianGradedArray` and drop the redundant per-block `haskey` guard, since matching axes imply matching allowed-block keys. Replace the `+`-only `mapreduce` method with a direct `sum`, the reduction actually needed, which avoids the footgun of skipping forbidden (zero) blocks for reductions where `f(0)` is nonzero.
Remove the unused `view(::AbstractGradedArray{T,1}, ::Block{1})` method, which rebuilt the same `Block{1}` and re-dispatched. It is shadowed by the concrete per-type `Block{1}` views and would otherwise recurse. Rename `axes_match_for_contraction` to `are_axes_contractible`.
Fix the `FusedGradedVector` docstring to say the stored blocks match the axis sectors exactly, which the constructor enforces. Add `dot` and `sum` tests, trim the unmatricize regression comment, and note that block-aware matmul and adjoint on `AbelianGradedMatrix` are still needed.
Route the block-wise `+`, `-`, `*`, and `/` on `FusedGradedMatrix` through a single `_broadcast_fusedgradedmatrix` helper, renamed from `_block_combine`. The binary form for `+` and `-` checks axis agreement and combines blocks over the union of sectors, and the single-array form applies a function to each stored block for the scalar ops. The name and a comment make clear that these methods stand in for the broadcasting currently disabled on `FusedGradedMatrix`, and are superseded once structure-preserving broadcasting lands.
Contraction now requires every contracted axis pair to be a canonical dual pair (`dual(ax1) == ax2`) for all sector types, fermionic included, dropping the earlier relaxation that also accepted same-`isdual` fermionic pairings. A contraction always pairs a space with its dual. Two fermionic contraction tests had inadvertently relied on same-`isdual` pairings. They now test the fermionic twist through proper canonical contractions, where the two odd contracted indices appear in swapped order so a single transposition introduces the expected -1.
TensorAlgebra 0.9.5 is registered with the projectto! and trivialrange primitives this PR overloads, so resolve it from the registry and set it as the compat floor.
| # axis types (`SectorRange`/`GradedOneTo`/`SectorOneTo`). Without the axis flip, | ||
| # `conj(t)` would leave bra-layer tensors with the same duality as the ket, and | ||
| # any contraction between them would silently pair non-dual against non-dual. | ||
| function Base.conj(a::AbelianGradedArray{T, N, D, S}) where {T, N, D, S} |
There was a problem hiding this comment.
This function does break the fermions again - it should be equal to permutedims(a', reverse(1:ndims(a)) which introduces -1 in the permutations
There was a problem hiding this comment.
I see, this one I definitely wasn't sure about, I should have run that by you before merging, happy to discuss and make a followup PR if needed.
| # Block-aware `LinearAlgebra.isdiag` for graded matrices. The generic | ||
| # `LinearAlgebra.isdiag` would fall through to `_isbanded_impl`'s scalar-indexed | ||
| # `iszero(view)` iteration, which throws on block storage. A graded matrix is | ||
| # diagonal iff it is block-diagonal (no off-diagonal blocks stored) and each | ||
| # stored block is itself diagonal — the latter checked block-by-block to avoid | ||
| # materializing the whole matrix. | ||
|
|
There was a problem hiding this comment.
Yes, good catch, thanks.
| twist(c::SectorRange) = TKS.twist(c.label) | ||
|
|
||
| Base.adjoint(r1::SectorRange) = dual(r1) | ||
| Base.conj(r1::SectorRange) = dual(r1) |
There was a problem hiding this comment.
Just wanted to double-check, what was the reason for changing this? (Also, this should probably be breaking, as this definitely breaks all scripts that were creating graded arrays, see the tests?)
For some background, there definitely are both the conjugate and the adjoint/dual representation, but there is a subtle difference, as the conjugate representation acts on the conjugate vector space, related to time-reversal things etc, by adjoint for the bra's, and not conj.
Obviously, we should not decide based solely on this, since realistically we never have a case where they don't match, but I also am not entirely sure about what this change really gains us, since the s' notation is a bit more compact and works just as well?
There was a problem hiding this comment.
The part that bothered me about using adjoint for this is that r' for a range r in Julia would return a column vector, i.e.:
julia> (1:3)'
1×3 adjoint(::UnitRange{Int64}) with eltype Int64:
1 2 3
i.e. as you know adjoint in Julia is both a conj and a transpose, so I was taking a (potentially simple-minded) view that we want adjoint without transpose, and therefore we want conj...
I didn't know all of that math background, indeed I wasn't totally sure if these definitions always coincided. I was taking a slightly more pedestrian view that I believed for our purposes they do coincide and for practical purposes in downstream packages like ITensorNetworksNext.jl it is a bit nicer to use conj instead of dual since conj is defined in Base (at least for now, my goal is for ITensorNetworksNext.jl to not depend on GradedArrays.jl, i.e. it is meant to be backend-agnostic. Obviously we could define dual in some shared interface but I didn't want to go there if possible).
I agree it is technically breaking but since no downstream packages are using it I considered it more of a bug fix/internal change, obviously that is up for debate and when the package is really being used I'd be more careful about marking breaking changes (i.e. the goal of this PR is to get ITensor/ITensorNetworksNext.jl#114 working with GradedArrays.jl, so after these sets of PRs ITensorNetworksNext.jl will be a proper "consumer" of GradedArrays.jl so if changes to GradedArrays breaks ITensorNetworksNext.jl I would consider that to be a breaking change).
There was a problem hiding this comment.
I think what is kind of interesting is that for a dual vector space, the 1x3 vector actually makes sense as well, it is the vector space of linear maps on vectors, i.e. things that take a vector and spit out a number, i.e. a 1x3 vector.
We then kind of lazily represent this by using a dual flag instead of an actual transpose range, but in principle this is all mathematically consistent.
From a math point of view, it is really the conj one that is the odd one out in a sense, for example A + conj(A) as matrices is a slightly strange operation, since A : V -> W means conj(A) : conj(V) -> conj(W) so you are kind of adding spaces that don't match, while A + A' is fine for square inputs, since A : V -> W results in A' : W -> V so if V == W (square) then this is well-defined.
I do absolutely agree that this is all kind of pedantic and really not super useful to make that a thing, but I would say that conj(V) is less convenient than V', which was a good argument for me 😉
There was a problem hiding this comment.
I see, I guess I wasn't thinking about the dual flag as effectively including a transposition as well.
Summary
Base.conjon sectors and axes (replaces the priorBase.adjoint), and onAbelianGradedArray(conjugates the data, flips axis duality).AbelianGradedArray:Base.copy,Base.copyto!,Base.fill!,Base.iszero,FI.zero!,scale!,Random.rand!/Random.randn!and theBase.rand/Base.randngraded-axes constructors, scalar*//,sum,Base.Array,LinearAlgebra.norm,LinearAlgebra.dot,LinearAlgebra.normalize.FusedGradedMatrix: block-wise+and-, scalar*//,LinearAlgebra.dot,LinearAlgebra.isdiag, and base matrix functions (sqrt,exp,log, ...). Elementwise broadcasting onFusedGradedMatrixis disabled, since the block structure makes it ill-defined. OnFusedGradedVector:Base.mapreduce,Base.map,MatrixAlgebraKit.diagonal.MatrixAlgebraKit.truncate(svd_trunc!, ::NTuple{3, FusedGradedMatrix}, ...)drops fully-truncated coupled sectors from the bond axis.TensorAlgebra.projectto!andTensorAlgebra.checked_projectto!onAbelianGradedArray, plusBlock{1}and graded-axisBase.similardisambiguators.TensorAlgebra.trivialrange(::Type{<:GradedOneTo}), the identity range fortensor_product_axison graded ranges (charge-0 one-dimensional sector).FI.permuteddimsto an eagerpermutedimsfallback at the abstract-array level (replacing the per-type methods), and add scalar (rank-0)unmatricize.GradedOneTohandling inmortar_axis, the matrix functions onFusedGradedMatrix, and==/hash.check_input(::typeof(contract), ::AbstractGradedArray, ...)requires every contracted-axis pair to be a canonical dual pair (dual(ax1) == ax2) for all sector types, throwing on same-isdualpairings that would otherwise silently align blocks index-wise.