[BlockSparseArrays] BlockSparseArray functionality

[ogauthe]

This issue lists functionalities and feature requests for BlockSparseArray.

using LinearAlgebra
using NDTensors.BlockSparseArrays: BlockSparseArray

a = BlockSparseArray{Float64}([2, 3], [2, 3]);

Bugs

• LinearAlgebra.Adjoint{T,<:BlockSparseArray{T}} returns a BlockedArray in certain slicing operations ([BlockSparseArrays] BlockSparseArray functionality #2).
• Fix cat(::BlockSparseArray, ::BlockSparseArray) for dual axes (followup to [BlockSparseArrays] Direct sum/cat ITensors.jl#1579).

Feature requests

• Allow the syntax BlockSparseArray{Float64}([U1(0) => 2, U1(1) => 3], [U1(0) => 2, U1(1) => 3]) which could implicitly create axes with GradedUnitRange internally.
• Use NestedPermutedDimsArray instead of SparsePermutedDimsArrayBlocks (similar to how we are removing SparseAdjointBlocks/SparseTransposeBlocks in [BlockSparseArrays] Simplifications of blocks for blocksparse Adjoint and Transpose ITensors.jl#1580). Started in [NDTensors] Introduce NestedPermutedDimsArrays submodule ITensors.jl#1589, [SparseArrayInterface] NestedPermutedDimsArray support ITensors.jl#1590.
• An alternative design to [NDTensors] Introduce NestedPermutedDimsArrays submodule ITensors.jl#1589 would be to redefine NestedPermutedDimsArray as a PermutedDimsArray wrapping a MappedArray where the map and inverse map convert to PermutedDimsArray, that would be good to explore so we don't have to support all of the NestedPermutedDimsArrays code, which is mostly just a copy of Base.PermutedDimsArrays anyway.
• If slices are just blockwise, like b = @view a[Block.(1:2), [Block(2), Block(1)], define blocks(b) as @view blocks(a)[1:2, [2, 1]], as opposed to using the more general SparseSubArrayBlocks in those cases. Like the new NestedPermutedDimsArray, in principle SparseSubArrayBlocks could be replaced by a NestedSubArray type that defines the slicing behavior of the array storing the blocks and also the slicing of the blocks themselves, but that might be overkill and the concept is very particular to block arrays. But maybe SubArray of the blocks could still be used to simplify the code logic in SparseSubArrayBlocks.
• Constructor from Dictionary should check block sizes ([BlockSparseArrays] BlockSparseArray functionality #2).
• Blockwise linear algebra operations like svd, qr, etc. See [ENHANCEMENT] Blockwise matrix factorizations #3. These are well defined if the block sparse matrix has a block structure (i.e. the sparsity pattern of the sparse array of arrays blocks(a)) corresponding to a generalized permutation matrix. Probably they should be called something like block_svd, block_eigen, block_qr, etc. to distinguish that they are meant to be used on block sparse matrices with those structures (and error if they don't have that structure). See ¹ for a prototype of a blockwise QR. See also BlockDiagonals.jl for an example in Julia of blockwise factorizations, they use a naming scheme svd_blockwise, eigen_blockwise, etc. The slicing operation introduced in [BlockSparseArrays] Sub-slices of multiple blocks ITensors.jl#1489 will be useful for performing block-wise truncated factorizations.
• Define storedblockview(a, ::Block) to get the view of a stored block, which would allow avoiding wrapping the block in a view. Also define getstoredindex(a, ::Block).
• Change the behavior of slicing with non-blocked ranges (such as a[1:2, 1:2]) to output non-blocked arrays, and define @blocked a[1:2, 1:2] to explicitly preserve blocking. See the discussion in Functionality for slicing with unit ranges that preserves block information JuliaArrays/BlockArrays.jl#347.
• Reconsider the design of how duality is stored in graded unit ranges (graded axes), for example storing it at the level of the sector labels with a new SectorDual type and/or as a boolean flag.
• Display non-initialized blocks differently from zero-initialized blocks, currently they both print as zeros.

Fixed

• Rename SparseArrayInterface to SparseArraysBase and move SparseArrayDOKs into the new SparseArraysBase. (Fixed in [SparseArraysBase] Rename SparseArrayInterface to SparseArraysBase ITensors.jl#1591, [SparseArraysBase] Absorb SparseArrayDOKs ITensors.jl#1592.)
• Constructors like BlockSparseArray{Float64,2,Matrix{Float64}}([2, 3], [2, 3]), BlockSparseArray{Float64,2}([2, 3], [2, 3]), BlockSparseMatrix{Float64}([2, 3], [2, 3]) are not defined. (Fixed by [BlockSparseArrays] Define more constructors ITensors.jl#1586.)
• Rename block_nstored to block_stored_length and nstored to stored_length. (Fixed by [SparseArrayInterface] [BlockSparseArrays] Rename nstored to stored_length ITensors.jl#1585.)
• Rename BlockSparseArrayLike to AnyAbstractBlockSparseArray, which is the naming convention used in other Julia packages for a similar concept².
• Implement direct sums/concatenations of block sparse arrays that preserve and make use of block sparsity. This could be done by overloading Base.cat and related functions. (Implemented in [BlockSparseArrays] Direct sum/cat ITensors.jl#1579.)
• Fix issues with nested slices like @views a[[Block(2), Block(1)], [Block(2), Block(1)]][2:4, 2:4] in Julia 1.11 (see tests marked as broken in [NDTensors] [ITensors] Update tests to use Julia version 1.10 and 1.11 ITensors.jl#1539). (Fixed in [BlockSparseArrays] Fix nested slicing in Julia 1.11 ITensors.jl#1575.)
• Support for blocks that are on GPU. (Partially addressed in [BlockSparseArrays] Initial support for more general blocks, such as GPU blocks ITensors.jl#1560 but more work is needed, we can track issues individually from now on.)
• Better support for blocks that are not Array, for example DiagonalArrays.DiagonalArray, SparseArrayDOKs.SparseArrayDOK, LinearAlgebra.Diagonal, etc. BlockSparseArray can have blocks that are AbstractArray subtypes, however some operations don't preserve those types properly (i.e. implicitly convert to Array blocks) or don't work. (Partially addressed in [BlockSparseArrays] Initial support for more general blocks, such as GPU blocks ITensors.jl#1560 but more work is needed, we can track issues individually from now on.)
• Cannot create zero dimension array (BlockSparseArray{Float64}() fails). (Fixed in [BlockSparseArrays] Zero dimensional block sparse array and some fixes for Adjoint and PermutedDimsArray ITensors.jl#1574.)
• permutedims crashes for some block sparse arrays ([BlockSparseArrays] BlockSparseArray functionality #2). (Fixed in [BlockSparseArrays] Zero dimensional block sparse array and some fixes for Adjoint and PermutedDimsArray ITensors.jl#1574.)
• copy(adjoint) does not preserve dual axes. (Fixed in [BlockSparseArrays] Zero dimensional block sparse array and some fixes for Adjoint and PermutedDimsArray ITensors.jl#1574.)
• block_stored_indices(::LinearAlgebra.Adjoint{T, BlockSparseArray}) does not transpose its indices. (Matt: I don't see an issue here, the values of block_stored_indices(a') are the nonzero/stored block locations, the keys of block_stored_indices(a') are an implementation detail and should not be used.)
• LinearAlgebra.norm(a) crashes when a contains NaN.
• Slicing a BlockSparseArray that has GradedUnitRange axes (as opposed to GradedOneTo) fails.
• Sub-slices of multiple blocks sometimes fails with a dual axis..
• a[:, :] creates an array with ill behaved axes (it should just be equivalent to copy(a)). Also triggers display error.
• Printing a BlockSparseArray triggers BoundsError in some rare contexts ([BlockSparseArrays] BlockSparseArray functionality #2).
• Automatically initializing a block when only a slice is used (this should work if you make sure to use a view, see [BlockSparseArrays] BlockSparseArray functionality #2).
• Implement slicing syntax for merging blocks, i.e. More convenient syntax for merging blocks JuliaArrays/BlockArrays.jl#359.
• a = BlockSparseArray{Float64}([2, 3], [2, 3]); @view a[Block(1, 1)] returns a SubArray where the last type parameter which marks whether or not the slice supports faster linear indexing is false, while it should be true if that is the case for that block of a (this is addressed by [BlockSparseArrays] Redesign block views again ITensors.jl#1513, @view a[Block(1, 1)] no longer outputs a SubArray, but rather either the block data directly or a BlockView object if the block doesn't exist yet).
• TensorAlgebra.contract fails when called with view(::BlockSparseArray, ::Block) or reshape(view(::BlockSparseArray, ::Block), ...). As a workaround it will work if you use view!/@view! instroduced in [BlockSparseArrays] Define in-place view that may instantiate blocks ITensors.jl#1498.
• Fix tests that are marked as broken in https://github.com/ITensor/ITensors.jl/blob/main/NDTensors/src/lib/BlockSparseArrays/test/test_basics.jl.
• a = BlockSparseArray{Float64}([2, 3], [2, 3]); b = @view a[Block.(1:2), Block.(1:2)]; b[Block(1, 1)] = randn(2, 2) doesn't set the block Block(1, 1) (it remains uninitialized, i.e. structurally zero). I think the issue is that @view b[Block(1, 1)] makes two layers of SubArray wrappers instead of flattening down to a single layer, and those two layers are not being dispatched on properly (in general we only catch if something is a BlockSparseArray or a BlockSparseArray wrapped in a single wrapper layer).
• Compatibility issues with BlockArrays v1.1, see CI for [SymmetrySectors] Non-abelian fusion ITensors.jl#1363. Fixed by [BlockSparseArrays] Update to BlockArrays v1.1, fix some issues with nested views ITensors.jl#1503.
• Rewrite GradedAxes using BlockArrays v1³.
• r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(r, r); size(view(a, Block(1,1))[1:1,1:1]) returns a tuple of LabelledInteger instead of Int (see discussion, keep it that way at least for now).
• r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(dual(r), r); @view(a[Block(1, 1)])[1:1, 1:1] and other combinations of dual lead to method ambiguity errors.
• Implement slicing syntax for slicing multiple subsets of blocks, i.e. Indexing with Vector{<:BlockIndexRange{1}} JuliaArrays/BlockArrays.jl#358.
• dual is not preserved when adding/subtracting BlockSparseArrays, i.e. g = gradedrange([U1(0) => 1]); m = BlockSparseArray{Float64}(dual(g), g); isdual(axes(m + m, 1)) should be true but is false.
• Error printing views of blocks of BlockSparseArray with GradedUnitRange axes, i.e.r = gradedrange([U1(0) => 1]); a = BlockSparseArray{Float64}(r, r); @view a[Block(1, 1)].
• Change the design of slicing with unit ranges, i.e. a[2:4, 2:4], by using BlockArrays.BlockSlice.
• Fix a[Block(2), Block(2)] = randn(3, 3).
• Fix a[Block(2, 2)] .= 1.
• Fix @view(a[Block(1, 1)])[1:1, 1:1] = 1.
• Throw error if a[Block(1, 1)] = b if size(a[Block(1, 1)]) != size(b).
• Fix matrix multiplication of BlockSparseMatrix involving dual axes.
• Fix matrix multiplication involving adjointed BlockSparseMatrix, i.e. a' * a and a * a', with and without dual axes.
• Take the dual of the axes in adjoint(::BlockSparseMatrix). Can be implemented by overloading axes(::Adjoint{<:Any,<:AbstractBlockSparseMatrix}).
• show(::Adjoint{<:Any,<:BlockSparseMatrix}) and show(::Transpose{<:Any,<:BlockSparseMatrix)) are broken.
• fix eachindex(::BlockSparseArray) involving dual axes.
• Fix [Hermitian] transposed BlockSparseMatrix, i.e. a' (in progress in⁴).
• Base.similar(a::BlockSparseArray, eltype::type) and Base.similar(a::BlockSparseArray, eltype::type, size::NTuple{N,AbstractUnitRange}) do not set eltype
• Make copy(::BlockSparseArray) copy the blocks.
• Slicing with a[1:2, 1:2] is not implemented yet and needs to be implemented (in progress in⁵).
• Function for accessing a list of initialized blocks. Currently you can use stored_indices(blocks(a))) to get a list of Block corresponding to initialized/stored blocks. Ideally there would be shorthands for this like block_stored_indices(a) (in progress in⁵).
• Function for accessing the number of initialized blocks. Currently you can use nstored(blocks(a)) to get the number of initialized/stored blocks. Ideally there would be shorthands for this like block_nstored(a) (in progress in⁵).
• In place operations .*= and ./=, such asa .*= 2, are broken (in progress in[^1]).
• Base.:*(::BlockSparseArray, x::Number) and Base.:/(::BlockSparseArray, x::Number) are not defined
• Base.:*(::ComplexF64, ::BlockSparseArray{Float64}) does not change data type for empty array and crashes if a contains data.

Footnotes

1. https://github.com/ITensor/ITensors.jl/blob/v0.3.57/NDTensors/src/lib/BlockSparseArrays/src/backup/LinearAlgebraExt/qr.jl ↩

2. https://github.com/JuliaGPU/GPUArrays.jl/blob/v11.1.0/lib/GPUArraysCore/src/GPUArraysCore.jl#L27, https://github.com/JuliaGPU/CUDA.jl/blob/v5.4.2/src/array.jl#L396 ↩

3. https://github.com/ITensor/ITensors.jl/pull/1452, https://github.com/JuliaArrays/BlockArrays.jl/pull/255 ↩

4. [BlockSparseArrays] Fix adjoint and transpose ITensors.jl#1470 ↩

5. [BlockSparseArrays] More general broadcasting and slicing ITensors.jl#1332 ↩ ↩² ↩³

[mtfishman]

Thanks. This currently works:

julia> using NDTensors.BlockSparseArrays: Block, BlockSparseArray, blocks

julia> using LinearAlgebra: I

julia> a = BlockSparseArray{Float64}([2, 2], [2, 2])
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
 0.0  0.0  │  0.0  0.0
 0.0  0.0  │  0.0  0.0
 ──────────┼──────────
 0.0  0.0  │  0.0  0.0
 0.0  0.0  │  0.0  0.0

julia> a[Block(2, 2)] = I(3)
3×3 Diagonal{Bool, Vector{Bool}}:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> a
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
 0.0  0.0  │  0.0  0.0
 0.0  0.0  │  0.0  0.0
 ──────────┼──────────
 0.0  0.0  │  1.0  0.0
 0.0  0.0  │  0.0  1.0

julia> using NDTensors.SparseArrayInterface: stored_indices

julia> stored_indices(blocks(a))
1-element Dictionaries.MappedDictionary{CartesianIndex{2}, CartesianIndex{2}, NDTensors.SparseArrayInterface.var"#1#2"{NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}}, Tuple{Dictionaries.Indices{CartesianIndex{2}}}}
 CartesianIndex(2, 2) │ CartesianIndex(2, 2)

though using this alternative syntax is currently broken:

julia> a = BlockSparseArray{Float64}([2, 2], [2, 2])
2×2-blocked 4×4 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{Vector{Int64}}, BlockArrays.BlockedUnitRange{Vector{Int64}}}}:
 0.0  0.0  │  0.0  0.0
 0.0  0.0  │  0.0  0.0
 ──────────┼──────────
 0.0  0.0  │  0.0  0.0
 0.0  0.0  │  0.0  0.0

julia> a[Block(2), Block(2)] = I(3)
ERROR: DimensionMismatch: tried to assign (3, 3) array to (2, 2) block
Stacktrace:
 [1] setindex!(::BlockSparseArray{…}, ::Diagonal{…}, ::Block{…}, ::Block{…})
   @ BlockArrays ~/.julia/packages/BlockArrays/L5yjb/src/abstractblockarray.jl:165
 [2] top-level scope
   @ REPL[30]:1
Some type information was truncated. Use `show(err)` to see complete types.

I would have to think about if it makes sense to support a[Block(2, 2)] .= 1.0 or a[Block(2), Block(2)] .= 1.0 if Block(2, 2) isn't initialized yet (probably we should support that).

[mtfishman]

In terms of svd and qr of BlockSparseMatrix, it is well defined to perform them block-wise (i.e. map the problem to factorizing the individual non-zero blocks) as long as there is only one block per row or column, i.e. diagonal up to block permutations of the rows and columns. So I was thinking we would have specialized block_svd and block_qr that checks if they have that structure and performs svd or qr block-wise, and otherwise it errors (it can check if blocks(a) is a generalized permutation matrix).

I have a protype of a QR decomposition of a BlockSparseMatrix here: https://github.com/ITensor/ITensors.jl/blob/v0.3.57/NDTensors/src/lib/BlockSparseArrays/src/backup/LinearAlgebraExt/qr.jl but it may be outdated since it was written based on an earlier version of BlockSparseArrays.

[mtfishman]

Also note that slicing like this should work right now:

a[Block(1, 1)[1:2, 1:2]]

i.e. you can take slices within a specified block. See BlockArrays.jl for a reference on that slicing notation.

[ogauthe]

new feature request: Base.:*(::BlockSparseArray, x::Number) and Base.:/(::BlockSparseArray, x::Number)

I updated the first comment.

Edit: FIXED

[ogauthe]

new issue: 1im * a does not modify a data type if a is empty. It crashes if a contains data.

Edit: FIXED

[ogauthe]

new issue: similar(a::BlockSparseArray, eltype::type) do not set new eltype.
similar(a::BlockSparseArray, eltype::type, size) behavior depends on size type.

• For size::NTuple{N,Int}, it returns a dense julia Array with correct eltype.
• For size::NTuple{N,AbstractUnitRange}, it returns a BlockSparseArray{Float64} and eltype is not correct.

edit: FIXED

[mtfishman]

@ogauthe a number of these issues were fixed by ITensor/ITensors.jl#1332, I've updated the list in the first post accordingly. I added regression tests in ITensor/ITensors.jl#1360 for ones that still need to be fixed, and additionally added placeholder tests that I've marked as broken in the BlockSparseArrays tests. Please continue to update this post with new issues you find, and/or make PRs with broken behavior marked with @test_broken or bug fixes and regression tests. I think with the reorganization I did in ITensor/ITensors.jl#1332 these kinds of issues will be easier to fix going forward, though I have to admit it may be a bit tough to jump into that code right now since there are many abstraction layers involved, since the goal is for that code to be as minimal and general as possible.

[ogauthe]

Feature request: display(::Adjoint). The function Base.adjoint is well defined and returns an Adjoint type, but it throws an error at display.
Closely related, it is not clear to me how GradedAxes.dual and BlockSparseArray should work together. The current behavior is to only accept a GradedUnitRange as an axis and to throw for a UnitRangeDual input. Base.adjoint does not dual the axis. I do not have an opinion yet whether this should be modified or not.

Edit: FIXED

[mtfishman]

I think ideally BlockSparseArray would accept any AbstractUnitRange.

Alternatively, BlockArrays.jl introduced an AbstractBlockedUnitRange in JuliaArrays/BlockArrays.jl#348, we could define a BlockedUnitRangeDual <: AbstractBlockedUnitRange type which wraps another AbstractBlockedUnitRange. I think it will require some experimentation to see which one leads to simpler code.

Good question about whether or not the axes should get dualed if Base.adjoint(::BlockSparseMatrix) is called. Probably we should dual the axes. For example, if we want operations like a' * a to work for a::BlockSparseMatrix and GradedUnitRange axes, I think we need to dual the axes.

[ogauthe]

The solution to accept any AbstractUnitRange is to replace the definition of Axes with

  Axes<:Tuple{Vararg{<:AbstractUnitRange,N}},

Then one can construct a BlockSparseArray with a dual axis. For a reason I do not really understand, the resulting array can be printed but not displayed. The error is not the same as when one tries to display a transpose or adoint.

g1 = gradedrange([U1(0) => 1])
m1 = BlockSparseArray{Float64}(dual(g1), g1,)

outputs

1×1-blocked 1×1 BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{UnitRangeDual{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}}}, Tuple{UnitRangeDual{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}}:
Error showing value of type BlockSparseArray{Float64, 2, Matrix{Float64}, NDTensors.SparseArrayDOKs.SparseArrayDOK{Matrix{Float64}, 2, NDTensors.BlockSparseArrays.BlockZero{Tuple{UnitRangeDual{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}}}, Tuple{UnitRangeDual{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}, BlockedUnitRange{Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}}}}:
ERROR: MethodError: no method matching NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}(::Int64)

Closest candidates are:
  (::Type{NDTensors.LabelledNumbers.LabelledInteger{Value, Label}} where {Value<:Integer, Label})(::Any, ::Any)
   @ NDTensors ~/Documents/itensor/ITensors.jl/NDTensors/src/lib/LabelledNumbers/src/labelledinteger.jl:2
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::BigFloat) where T<:Integer
   @ Base mpfr.jl:378
  ...

Stacktrace:
  [1] convert(::Type{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}, x::Int64)
    @ Base ./number.jl:7
  [2] cvt1
    @ ./essentials.jl:468 [inlined]
  [3] ntuple
    @ ./ntuple.jl:49 [inlined]
  [4] convert(::Type{Tuple{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}, Int64}}, x::Tuple{Int64, Int64})
    @ Base ./essentials.jl:470
  [5] push!(a::Vector{Tuple{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}, Int64}}, item::Tuple{Int64, Int64})
    @ Base ./array.jl:1118
  [6] alignment(io::IOContext{Base.TTY}, X::AbstractVecOrMat, rows::Vector{NDTensors.LabelledNumbers.LabelledInteger{Int64, U1{Int64}}}, cols::Vector{Int64}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Int64)
    @ Base ./arrayshow.jl:76
  [7] _print_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::UnitRange{…}, colsA::UnitRange{…})
    @ Base ./arrayshow.jl:207
  [8] print_matrix(io::IOContext{…}, X::BlockSparseArray{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base ./arrayshow.jl:171
  [9] print_matrix
    @ ./arrayshow.jl:171 [inlined]
 [10] print_array
    @ ./arrayshow.jl:358 [inlined]
 [11] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::BlockSparseArray{Float64, 2, Matrix{…}, NDTensors.SparseArrayDOKs.SparseArrayDOK{…}, Tuple{…}})
    @ Base ./arrayshow.jl:399
 [12] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:273
 [13] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [14] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:259
 [15] display
    @ ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:278 [inlined]
 [16] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:340
 [17] #invokelatest#2
    @ ./essentials.jl:892 [inlined]
 [18] invokelatest
    @ ./essentials.jl:889 [inlined]
 [19] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:315
 [20] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:284
 [21] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [22] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:282
 [23] (::REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:911
 [24] #invokelatest#2
    @ ./essentials.jl:892 [inlined]
 [25] invokelatest
    @ ./essentials.jl:889 [inlined]
 [26] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/LineEdit.jl:2656
 [27] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:1312
 [28] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ~/.julia/juliaup/julia-1.10.3+0.x64.linux.gnu/share/julia/stdlib/v1.10/REPL/src/REPL.jl:386
Some type information was truncated. Use `show(err)` to see complete types.

Edit: FIXED

[mtfishman]

Thanks for investigating. That seems like the right move to generalize the axes in that way. Hopefully that error is easy enough to circumvent.

[ogauthe]

I continue in exploring the effect of dual. eachindex is broken on such an array. This affects ==.

julia> eachindex(m1)
ERROR: MethodError: no method matching AbstractUnitRange{…}(::NDTensors.GradedAxes.UnitRangeDual{…})

Closest candidates are:
  AbstractUnitRange{T}(::AbstractUnitRange{T}) where T
   @ Base range.jl:1308
  AbstractUnitRange{T}(::NDTensors.LabelledNumbers.LabelledUnitRange) where T
   @ NDTensors ~/Documents/itensor/ITensors.jl/NDTensors/src/lib/LabelledNumbers/src/labelledunitrange.jl:17
  AbstractUnitRange{Int64}(::Static.OptionallyStaticUnitRange)
   @ Static ~/.julia/packages/Static/6JeQI/src/ranges.jl:258
  ...

Stacktrace:
 [1] OrdinalRange{…}(r::NDTensors.GradedAxes.UnitRangeDual{…})
   @ Base ./range.jl:1313
 [2] convert(::Type{…}, r::NDTensors.GradedAxes.UnitRangeDual{…})
   @ Base ./range.jl:265
 [3] (::Base.IteratorsMD.var"#3#4")(r::NDTensors.GradedAxes.UnitRangeDual{…})
   @ Base.IteratorsMD ./multidimensional.jl:259
 [4] map
   @ ./tuple.jl:292 [inlined]
 [5] CartesianIndices(inds::Tuple{NDTensors.GradedAxes.UnitRangeDual{…}, BlockArrays.BlockedUnitRange{…}})
   @ Base.IteratorsMD ./multidimensional.jl:259
 [6] eachindex(::IndexCartesian, A::BlockSparseArray{…})
   @ Base.IteratorsMD ./multidimensional.jl:392
 [7] eachindex(A::BlockSparseArray{Float64, 2, Matrix{…}, NDTensors.SparseArrayDOKs.SparseArrayDOK{…}, Tuple{…}})
   @ Base ./abstractarray.jl:378
 [8] top-level scope
   @ REPL[37]:1
Some type information was truncated. Use `show(err)` to see complete types.

Edit: FIXED