From 44b97de6f5e9df31f5a822b2212ec7d60e6ab72e Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Wed, 8 Jul 2026 13:45:27 -0400 Subject: [PATCH 1/2] Bare-matrix factorizations, *, and one! for AbelianGradedMatrix Route the matrix-level linear algebra for a bare `AbelianGradedMatrix` through the matricizing `TensorAlgebra` factorizations: the `MatrixAlgebraKit` factorization family, `Base.:*`, and in-place `one!`. --- Project.toml | 6 +- src/abeliangradedarray.jl | 12 ++++ src/matrixalgebrakit.jl | 24 ++++++++ test/test_factorizations.jl | 112 +++++++++++++++++++++++++++++++++++- test/test_tensoralgebra.jl | 14 ++--- 5 files changed, 157 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index a99bd82e..31e3208a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "GradedArrays" uuid = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2" -version = "0.13.8" +version = "0.13.9" authors = ["ITensor developers and contributors"] [workspace] @@ -26,6 +26,10 @@ TensorKitSectors = "13a9c161-d5da-41f0-bcbd-e1a08ae0647f" SUNRepresentations = "1a50b95c-7aac-476d-a9ce-2bfc675fc617" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" +[sources.TensorAlgebra] +rev = "mf/nextgen-followups-round1" +url = "https://github.com/ITensor/TensorAlgebra.jl" + [extensions] GradedArraysSUNRepresentationsExt = "SUNRepresentations" GradedArraysTensorKitExt = "TensorKit" diff --git a/src/abeliangradedarray.jl b/src/abeliangradedarray.jl index d627f321..c9ee26a8 100644 --- a/src/abeliangradedarray.jl +++ b/src/abeliangradedarray.jl @@ -822,3 +822,15 @@ end function Base.getindex(a::Array, ax1::GradedOneTo) return invoke(getindex, Tuple{AbstractArray, GradedOneTo, Vararg{GradedOneTo}}, a, ax1) end + +# --------------------------------------------------------------------------- +# matrix multiplication +# --------------------------------------------------------------------------- + +# Matrix-matrix multiply via `TensorAlgebra.contract` (which matricizes internally). The +# contracted label `2` pairs `a`'s column axis with `b`'s row axis; the explicit output +# labels `(1, 3)` fix the result as `a`'s row axis (codomain) × `b`'s column axis (domain), +# so `(a * b)[i, j] == sum_k a[i, k] * b[k, j]`. +function Base.:*(a::AbelianGradedMatrix, b::AbelianGradedMatrix) + return TensorAlgebra.contract((1, 3), a, (1, 2), b, (2, 3)) +end diff --git a/src/matrixalgebrakit.jl b/src/matrixalgebrakit.jl index 0c4e03d1..b2fe25fb 100644 --- a/src/matrixalgebrakit.jl +++ b/src/matrixalgebrakit.jl @@ -34,6 +34,30 @@ for f in [ end end +# Bare-matrix factorizations +# --------------------------- +# There is no in-place block algorithm for an unfused `AbelianGradedMatrix`, so the plain +# matrix forms (`MAK.svd_compact(m)`, etc.) route through the matricizing `TensorAlgebra` +# factorizations: matricize to a `FusedGradedMatrix`, run the block factorization, then +# unmatricize back. The factors are returned as graded matrices. Only the functions with an +# identically-named `TensorAlgebra` perm-form are delegated here (`qr_null`/`lq_null` are +# spelled `left_null`/`right_null` there, and `project_antihermitian`/`project_isometric` +# have no perm-form). +for f in ( + :svd_compact, :svd_full, :svd_vals, :qr_compact, :qr_full, :lq_compact, + :lq_full, :eig_full, :eig_vals, :eigh_full, :eigh_vals, :left_polar, + :right_polar, :project_hermitian, + ) + @eval function MAK.$f(m::AbelianGradedMatrix; kwargs...) + return TensorAlgebra.$f(m, (1,), (2,); kwargs...) + end +end + +# In-place graded identity fill. Filling the unfused data blocks with identities is not the +# graded identity map in general, so route through the fused path: `TensorAlgebra.one!` +# matricizes, fills the fused matrix with `MAK.one!`, and scatters it back into `a`. +MAK.one!(a::AbelianGradedMatrix) = TensorAlgebra.one!(a, Val(1)) + # Generic Implementations # ----------------------- # utility function to do something with each block diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 2c89715a..1ea8e72b 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,7 +1,9 @@ import MatrixAlgebraKit as MAK -using GradedArrays: FusedGradedMatrix, FusedGradedVector, GradedBlockAlgorithm, U1, Z2 +using GradedArrays: AbelianGradedMatrix, FusedGradedMatrix, FusedGradedVector, + GradedBlockAlgorithm, U1, Z2, dual, gradedrange using LinearAlgebra: Diagonal, I, eigvals, isposdef, istril, istriu, norm using MatrixAlgebraKit: isisometric, isunitary +using TensorAlgebra: TensorAlgebra using Test: @test, @testset # --------------------------------------------------------------------------- @@ -408,4 +410,112 @@ end end end end + + # ----------------------------------------------------------------------- + @testset "matrix multiplication" begin + g = gradedrange([U1(0) => 2, U1(1) => 3, U1(2) => 2]) + h = gradedrange([U1(0) => 3, U1(1) => 2, U1(2) => 4]) + a = randn(Float64, (g, dual(g))) + b = randn(Float64, (g, dual(h))) + c = a * b + @test c isa AbelianGradedMatrix + # `(a * b)[i, j] == sum_k a[i, k] * b[k, j]`. + @test Array(c) ≈ Array(a) * Array(b) + # Result axes: codomain from `a`, domain from `b`. + @test axes(c, 1) == axes(a, 1) + @test axes(c, 2) == axes(b, 2) + end + + # ----------------------------------------------------------------------- + # Bare-matrix factorizations delegate to the matricizing `TensorAlgebra` forms. + @testset "factorizations on a bare AbelianGradedMatrix" begin + g = gradedrange([U1(0) => 2, U1(1) => 3, U1(2) => 2]) + h = gradedrange([U1(0) => 3, U1(1) => 2, U1(2) => 4]) + m_rect = randn(Float64, (g, dual(h))) + m_sq = randn(Float64, (g, dual(g))) + m_herm = MAK.project_hermitian(randn(Float64, (g, dual(g)))) + + # helper: compare two `FusedGradedVector`s block-by-block (the broadcasting `-` + # path they would otherwise take is not supported). + fgv_approx(x, y) = + keys(x.blocks) == keys(y.blocks) && + all(x.blocks[k] ≈ y.blocks[k] for k in keys(x.blocks)) + + @testset "svd_compact" begin + U, S, Vᴴ = MAK.svd_compact(m_rect) + @test all(x -> x isa AbelianGradedMatrix, (U, S, Vᴴ)) + @test axes(U, 1) == axes(m_rect, 1) + @test axes(Vᴴ, 2) == axes(m_rect, 2) + @test U * S * Vᴴ ≈ m_rect + @test Array(U) * Array(S) * Array(Vᴴ) ≈ Array(m_rect) + end + + @testset "svd_full" begin + U, S, Vᴴ = MAK.svd_full(m_rect) + @test all(x -> x isa AbelianGradedMatrix, (U, S, Vᴴ)) + @test Array(U) * Array(S) * Array(Vᴴ) ≈ Array(m_rect) + end + + @testset "svd_vals" begin + @test fgv_approx(MAK.svd_vals(m_rect), MAK.svd_vals(FusedGradedMatrix(m_rect))) + end + + @testset "qr_compact / qr_full" begin + Q, R = MAK.qr_compact(m_rect) + @test Array(Q) * Array(R) ≈ Array(m_rect) + Q, R = MAK.qr_full(m_rect) + @test Array(Q) * Array(R) ≈ Array(m_rect) + end + + @testset "lq_compact / lq_full" begin + L, Q = MAK.lq_compact(m_rect) + @test Array(L) * Array(Q) ≈ Array(m_rect) + L, Q = MAK.lq_full(m_rect) + @test Array(L) * Array(Q) ≈ Array(m_rect) + end + + @testset "eig_full / eig_vals" begin + D, V = MAK.eig_full(m_sq) + @test Array(m_sq) * Array(V) ≈ Array(V) * Array(D) + @test fgv_approx(MAK.eig_vals(m_sq), MAK.eig_vals(FusedGradedMatrix(m_sq))) + end + + @testset "eigh_full / eigh_vals" begin + D, V = MAK.eigh_full(m_herm) + @test Array(m_herm) ≈ Array(V) * Array(D) * Array(V)' + @test fgv_approx( + MAK.eigh_vals(m_herm), + MAK.eigh_vals(FusedGradedMatrix(m_herm)) + ) + end + + @testset "left_polar / right_polar" begin + W, P = MAK.left_polar(m_sq) + @test Array(W) * Array(P) ≈ Array(m_sq) + P, W = MAK.right_polar(m_sq) + @test Array(P) * Array(W) ≈ Array(m_sq) + end + + @testset "project_hermitian" begin + m = randn(Float64, (g, dual(g))) + @test Array(MAK.project_hermitian(m)) ≈ (Array(m) + Array(m)') / 2 + end + end + + # ----------------------------------------------------------------------- + @testset "one! on an AbelianGradedMatrix" begin + g = gradedrange([U1(0) => 2, U1(1) => 3, U1(2) => 2]) + a = randn(Float64, (g, dual(g))) + a_before = Array(copy(a)) + + b = MAK.one!(a) + # In place: returns `a` and mutates its contents. + @test b === a + @test Array(a) != a_before + + # Same as the existing graded identity constructor and the dense identity. + id = TensorAlgebra.one(randn(Float64, (g, dual(g))), (1,), (2,)) + @test Array(a) ≈ Array(id) + @test Array(a) ≈ Matrix(1.0I, size(a)...) + end end # @testset "Factorizations" diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl index c7ffdf07..d3763718 100644 --- a/test/test_tensoralgebra.jl +++ b/test/test_tensoralgebra.jl @@ -449,14 +449,11 @@ end A = AbelianGradedArray{Float64}(undef, s, dual(s)) randn!(A) U, S, Vᴴ = TensorAlgebra.svd_compact(A, (1,), (2,)) - # The natural `U * S * Vᴴ` form falls into LinearAlgebra's `_tri_matmul`, - # which scalar-indexes on `AbstractGradedArray`. `contract` is the - # block-wise route, but the chain form should also work once a block-aware - # matmul lands on `AbstractGradedMatrix`. US = contract((:a, :r), U, (:a, :i), S, (:i, :r)) USV = contract((:a, :b), US, (:a, :r), Vᴴ, (:r, :b)) @test A ≈ USV - @test_broken A ≈ U * S * Vᴴ + # `*` on `AbelianGradedMatrix` routes through the block-wise `contract`. + @test A ≈ U * S * Vᴴ end @testset "TA.gram_eigh_full_with_pinv on AbelianGradedMatrix (axes_Y regression)" begin @@ -467,14 +464,13 @@ end # so we stay on the graded matmul path; the natural `*` form is broken # against the same scalar-indexing path as the SVD round-trip above. A = contract((:a, :b), B, (:a, :r), conj(B), (:b, :r)) + # `*` on two `AbelianGradedMatrix` works, but the adjoint forms (`B * B'`, + # `X * X'` below) still need a block-aware `adjoint`; `B'` is an `Adjoint` + # wrapper that falls through to LinearAlgebra's scalar-indexing path. @test_broken A ≈ B * B' X, Y = TensorAlgebra.gram_eigh_full_with_pinv(A, (1,), (2,)) # X · conj(X) ≈ A on the rank subspace. @test A ≈ contract((:a, :b), X, (:a, :r), conj(X), (:b, :r)) - # Matmul (`*`) on `AbelianGradedMatrix` is unimplemented, so any `X * Y` - # falls through to LinearAlgebra's scalar-indexing path and throws. The - # adjoint forms (`X * X'`, `B * B'`) additionally need a block-aware - # `adjoint`. Both will pass once those land. @test_broken A ≈ X * X' # Y is a left inverse of X on the rank subspace. YX = contract((:r, :s), Y, (:r, :a), X, (:a, :s)) From f2c6b251b52075e1128e46fd3b34014db09b6ff7 Mon Sep 17 00:00:00 2001 From: Matthew Fishman Date: Thu, 9 Jul 2026 11:39:20 -0400 Subject: [PATCH 2/2] Add ungrade and matrix tr for graded and sector arrays Extends TensorAlgebra.ungrade for GradedOneTo, mapping a graded range to the plain range over its total dimension so a range and its dual share an ungraded value. Adds LinearAlgebra.tr for graded matrices and the sector block types. A graded matrix sums the traces of its diagonal blocks, and each block trace factorizes into the sector's quantum dimension times the trace of the reduced data. This is the plain matricized trace, with no fermionic twist sign. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/abstractgradedarray.jl | 10 ++++++++++ src/abstractsectorarray.jl | 9 +++++++++ src/abstractsectordelta.jl | 4 ++++ src/gradedoneto.jl | 4 ++++ test/test_gradedoneto.jl | 6 ++++++ test/test_sectoridentity.jl | 6 ++++++ test/test_sectormatrix.jl | 7 +++++++ test/test_tensoralgebra.jl | 13 +++++++++++++ 8 files changed, 59 insertions(+) diff --git a/src/abstractgradedarray.jl b/src/abstractgradedarray.jl index 0a4bdf3d..3e290c9c 100644 --- a/src/abstractgradedarray.jl +++ b/src/abstractgradedarray.jl @@ -15,6 +15,16 @@ function isblockdiagonal(A::AbstractGradedMatrix) return true end +# Trace: sum the traces of the diagonal blocks (same block position on both axes). Each block +# view carries its sector, so its `tr` includes that sector's quantum dimension. This is the +# plain matricized trace, with no fermionic twist sign. +function LinearAlgebra.tr(A::AbstractGradedMatrix) + return sum(eachblockstoredindex(A); init = zero(eltype(A))) do bI + row, col = Tuple(bI) + return row == col ? LinearAlgebra.tr(view(A, bI)) : zero(eltype(A)) + end +end + # Scalar indexing is not supported for graded arrays. function Base.getindex(::AbstractGradedArray, ::Vararg{Int}) return error( diff --git a/src/abstractsectorarray.jl b/src/abstractsectorarray.jl index 31cc24ce..5f337adf 100644 --- a/src/abstractsectorarray.jl +++ b/src/abstractsectorarray.jl @@ -71,6 +71,15 @@ function TensorAlgebra.zero!(a::AbstractSectorArray) return a end +# ======================== trace ======================== + +# A 2D sector block is the tensor product of its structural factor `sector(a)` and its reduced +# data `data(a)`, so the trace factorizes: the sector's quantum dimension times the trace of +# the reduced data. +function LinearAlgebra.tr(a::AbstractSectorArray{<:Any, <:Any, 2}) + return LinearAlgebra.tr(sector(a)) * LinearAlgebra.tr(data(a)) +end + # ======================== display ======================== function Base.print_array(io::IO, sa::AbstractSectorArray) diff --git a/src/abstractsectordelta.jl b/src/abstractsectordelta.jl index 7667736b..4208562c 100644 --- a/src/abstractsectordelta.jl +++ b/src/abstractsectordelta.jl @@ -13,3 +13,7 @@ sectortype(::Type{<:AbstractSectorDelta{T, S}}) where {T, S} = S Base.copy(A::AbstractSectorDelta) = A Base.size(A::AbstractSectorDelta) = length.(axes(A)) + +# A 2D structural delta on a diagonal (coupled) block is the identity over the sector, so its +# trace is the sector's quantum dimension: the length of the diagonal. +LinearAlgebra.tr(A::AbstractSectorDelta{<:Any, <:Any, 2}) = size(A, 1) diff --git a/src/gradedoneto.jl b/src/gradedoneto.jl index f2111dd9..29d86097 100644 --- a/src/gradedoneto.jl +++ b/src/gradedoneto.jl @@ -67,6 +67,10 @@ function TensorAlgebra.trivialrange(::Type{GradedOneTo{S}}, n::Integer) where {S return gradedrange([trivial(S) => n]) end +# The ungraded extent of a graded range is the plain range over its total dimension, dropping +# sectors and the arrow so a range and its `dual` share an ungraded value. +TensorAlgebra.ungrade(g::GradedOneTo) = Base.OneTo(length(g)) + """ gradedrange(xs::AbstractVector{<:Pair}) diff --git a/test/test_gradedoneto.jl b/test/test_gradedoneto.jl index de9318ba..75ce5949 100644 --- a/test/test_gradedoneto.jl +++ b/test/test_gradedoneto.jl @@ -230,4 +230,10 @@ using Test: @test, @test_throws, @testset g = gradedrange([U1(0) => 2, U1(1) => 3]) @test TensorAlgebra.to_range(g) === g end + + @testset "ungrade drops sectors and arrow" begin + g = gradedrange([U1(0) => 2, U1(1) => 3]) + @test TensorAlgebra.ungrade(g) == Base.OneTo(length(g)) + @test TensorAlgebra.ungrade(dual(g)) == TensorAlgebra.ungrade(g) + end end diff --git a/test/test_sectoridentity.jl b/test/test_sectoridentity.jl index 56ec239a..e3d83f1d 100644 --- a/test/test_sectoridentity.jl +++ b/test/test_sectoridentity.jl @@ -1,5 +1,6 @@ using GradedArrays: SU2, SectorIdentity, SectorRange, U1, dual, isdual, sectoraxes, sectortype +using LinearAlgebra: tr using TensorKitSectors: TensorKitSectors as TKS using Test: @test, @testset @@ -65,4 +66,9 @@ using Test: @test, @testset @test isdual(si, 1) == false @test isdual(si, 2) == true end + + @testset "tr — quantum dimension of the sector" begin + @test tr(SectorIdentity{Float64}(U1(0))) == 1 + @test tr(SectorIdentity{Float64}(SU2(1 // 2))) == 2 + end end diff --git a/test/test_sectormatrix.jl b/test/test_sectormatrix.jl index e7a0c29d..3d6b9856 100644 --- a/test/test_sectormatrix.jl +++ b/test/test_sectormatrix.jl @@ -1,6 +1,7 @@ using GradedArrays: GradedArrays, AbelianSectorArray, SU2, SectorIdentity, SectorMatrix, SectorOneTo, SectorRange, U1, data, dataaxes, dual, isdual, sector, sector_kron, sectoraxes, sectortype +using LinearAlgebra: tr using TensorKitSectors: TensorKitSectors as TKS using Test: @test, @test_throws, @testset @@ -143,4 +144,10 @@ using Test: @test, @test_throws, @testset @test sm isa SectorMatrix{Float64, U1, Matrix{Float64}} @test size(data(sm)) == (3, 4) end + + @testset "tr — sector quantum dimension times reduced-data trace" begin + d = [1.0 2.0; 3.0 4.0] + @test tr(SectorMatrix(U1(0), d)) == tr(d) # dim 1 + @test tr(SectorMatrix(SU2(1 // 2), d)) == 2 * tr(d) # dim 2 + end end diff --git a/test/test_tensoralgebra.jl b/test/test_tensoralgebra.jl index d3763718..852e769a 100644 --- a/test/test_tensoralgebra.jl +++ b/test/test_tensoralgebra.jl @@ -5,6 +5,7 @@ using GradedArrays: AbelianGradedArray, AbelianGradedMatrix, AbelianSectorArray, SectorRange, U1, data, datalengths, dual, eachblockstoredindex, eachsectoraxis, flip, gradedrange, isdual, sector, sectoraxes, sectormergesort, sectors, sectortype, tensor_product +using LinearAlgebra: tr using Random: randn! using TensorAlgebra: TensorAlgebra, FusionStyle, contract, linearbroadcasted, matricize, matricizeperm, unmatricize @@ -221,6 +222,18 @@ end @test data(fsm[Block(3, 3)]) ≈ 2 * ones(1, 1) end +@testset "tr of a matricized AbelianGradedArray" begin + g = gradedrange([U1(0) => 1, U1(1) => 1]) + a = ones(Float64, g, g, dual(g), dual(g)) + a[Block(2, 2, 2, 2)] .*= 2 # give the blocks distinct traces + + # `tr` on the matricized graded matrix sums the diagonal blocks and matches the dense trace. + fsm = matricizeperm(a, (1, 2), (3, 4)) + @test tr(fsm) ≈ tr(Array(fsm)) + # `TensorAlgebra.tr` over the (1, 2) | (3, 4) bipartition routes through the same path. + @test TensorAlgebra.tr(a, (1, 2, 3, 4), (1, 2), (3, 4)) ≈ tr(Array(fsm)) +end + @testset "matricize 3D AbelianGradedArray and unmatricize round-trip" begin # 3D case where the merged codomain (tensor product of two `r`s) has # sectors absent from the domain — the asymmetric `FusedGradedMatrix`