From 0f73d9db05567583a11fe31c23a88ec3664300b5 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Wed, 10 Jul 2024 14:58:42 +0100 Subject: [PATCH] Change sorting and sort young tableau (#7) * Change sorting and sort young tableau * sort young tableaux * Update runtests.jl * Take advantage of sparsity in Q --- Project.toml | 5 +- src/NumericalRepresentationTheory.jl | 138 ++++++++++++++++++++++----- test/runtests.jl | 78 ++++++++++----- 3 files changed, 175 insertions(+), 46 deletions(-) diff --git a/Project.toml b/Project.toml index c3f7e39..92c660f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,12 @@ name = "NumericalRepresentationTheory" uuid = "6b7c1d51-ecee-4149-97a8-50646b514dce" authors = ["Sheehan Olver "] -version = "0.2.0" +version = "0.3.0" [deps] +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Permutations = "2ae35dd2-176d-5d53-8349-f30d82d94d4f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/src/NumericalRepresentationTheory.jl b/src/NumericalRepresentationTheory.jl index 9c3e4af..122f76c 100644 --- a/src/NumericalRepresentationTheory.jl +++ b/src/NumericalRepresentationTheory.jl @@ -1,9 +1,10 @@ module NumericalRepresentationTheory -using Base, LinearAlgebra, Permutations, SparseArrays +using Base, LinearAlgebra, Permutations, SparseArrays, BlockBandedMatrices, BlockArrays, FillArrays import Base: getindex, size, setindex!, maximum, Int, length, ==, isless, copy, kron, hash, first, show, lastindex, |, Integer, BigInt import LinearAlgebra: adjoint, transpose, eigen +import SparseArrays: blockdiag ## Kronecker product of Sn @@ -40,13 +41,14 @@ Partition(σ::Int...) = Partition([σ...]) function isless(a::Partition, b::Partition) n,m = Int(a), Int(b) - n < m && return true - if n == m - M,N = length(a.σ), length(b.σ) - M ≠ N && return isless(M,N) - for k = 1:N - a.σ[k] > b.σ[k] && return true - a.σ[k] < b.σ[k] && return false + if n ≠ m + return n < m + end + + M,N = length(a.σ), length(b.σ) + for k = 1:min(M,N) + if a.σ[k] ≠ b.σ[k] + return a.σ[k] < b.σ[k] end end @@ -128,9 +130,9 @@ YoungMatrix(dat, σ::Vector{Int}) = YoungMatrix(dat, Partition(σ)) copy(Y::YoungMatrix) = YoungMatrix(copy(Y.data), Y.rows, Y.columns) size(Y::YoungMatrix) = size(Y.data) -getindex(Y::YoungMatrix, k::Int, j::Int) = ifelse(k ≤ Y.columns[j] && j ≤ Y.rows[k], Y.data[k,j], 0) +getindex(Y::YoungMatrix, k::Int, j::Int) = ifelse(k ≤ Y.columns[j] && j ≤ Y.rows[k], Y.data[k,j], 0) function setindex!(Y::YoungMatrix, v, k::Int, j::Int) - @assert k ≤ Y.columns[j] && j ≤ Y.rows[k] + @assert k ≤ Y.columns[j] && j ≤ Y.rows[k] Y.data[k,j] = v end @@ -141,6 +143,16 @@ struct YoungTableau partitions::Vector{Partition} end +function isless(a::YoungTableau, b::YoungTableau) + if length(a.partitions) ≠ length(b.partitions) + return length(a.partitions) < length(b.partitions) + end + for (p,q) in zip(a.partitions, b.partitions) + p ≠ q && return isless(p,q) + end + return false +end + function YoungMatrix(Yt::YoungTableau) ps = Yt.partitions @@ -189,7 +201,7 @@ end function youngtableaux(σ::Partition) σ == Partition([1]) && return [YoungTableau([σ])] Yts = mapreduce(youngtableaux, vcat, lowerpartitions(σ)) - map(Yt -> YoungTableau([Yt.partitions; σ]), Yts) + sort!(map(Yt -> YoungTableau([Yt.partitions; σ]), Yts)) end @@ -274,6 +286,14 @@ end Representation(σ::Int...) = Representation(Partition(σ...)) Representation(σ::Partition) = Representation(irrepgenerators(σ)) kron(A::Representation, B::Representation) = Representation(kron.(A.generators, B.generators)) + + +blockdiag(A::Representation...) = Representation(blockdiag.(getfield.(A, :generators)...)) + + +conjugate(ρ::Representation, Q) = Representation(map(g -> Q'*g*Q, ρ.generators)) + + ⊗(A::Representation, B::Representation) = kron(A, B) |(A::Representation, n::AbstractVector) = Representation(A.generators[n]) @@ -329,6 +349,23 @@ end gelfand_reduce(R::Representation) = gelfand_reduce(Matrix.(gelfandbasis(R.generators))) +function sortcontentsperm(Λ) + ps = contents2partition(Λ) + perm = Int[] + + + for pₖ in union(ps) + dₖ = hooklength(pₖ) + inds = findall(==(pₖ), ps) + aₖ = length(inds) ÷ dₖ # number of times repeated + for j = 1:aₖ + append!(perm, inds[j:aₖ:end]) + end + end + + perm +end + function gelfand_reduce(X) λ̃, Q₁ = eigen(Symmetric(X[1])) λ = round.(Int, λ̃) @@ -359,12 +396,23 @@ function gelfand_reduce(X) Λ, Q₁*Q end +""" + singlemultreduce(ρ) + +reduces a representation containing only a single irrep, multiple times. +""" + function singlemultreduce(ρ) m = multiplicities(ρ) @assert length(m) == 1 singlemultreduce(ρ, Representation(first(keys(m)))) end +""" + singlemultreduce(ρ, σ) + +reduces a representation containing only a single irrep `σ`, multiple times. +""" function singlemultreduce(ρ, σ) m = size(σ,1) n = size(ρ,1) @@ -372,11 +420,57 @@ function singlemultreduce(ρ, σ) Q̃ = nullspace(convert(Matrix,A); atol=1E-10)*sqrt(m) reshape(vec(Q̃), n, n) end + + +""" + singlemultreduce_blockdiag(ρ, σ) + +reduces a representation containing only a single irrep `σ`, multiple times. +Unlike `singlemultreduce`, it is assumed that `ρ` comes from `gelfand_reduce` and hence +the returned `Q` is guaranteed to be block-diagonal. +""" +function singlemultreduce_blockdiag(ρ, σ) + tol = 1E-10 + m = size(σ,1) + n = size(ρ,1) + μ = n ÷ m + A = vcat((kron.(Ref(I(m)), ρ.generators) .- kron.(σ.generators, Ref(I(n))))...) + + # determine columns of `A` corresponding to non-zero entries of `Q` + jr = Int[] + for k = 1:m + append!(jr, range((k-1)*μ*m+k; step=m, length=μ)) + end + + # determine non-zero rows of A[:,jr] + kr = Int[] + for k = 1:size(A,1) + if norm(A[k,jr]) > tol + push!(kr, k) + end + end + + V = nullspace(A[kr,jr])*sqrt(m) + + # now populate non-zero entries of `Q` + Q = BandedBlockBandedMatrix{Float64}(undef, Fill(m,μ), Fill(m,μ), (n-1,n-1), (0,0)) + for j = 1:size(V,2) + sh = (j-1) * size(V,1) + for k = 1:size(V,1) + view(Q,:,Block(j))[jr[k]] = V[k,j] + end + end + + Q +end function blockdiagonalize(ρ::Representation) Λ,Q = gelfand_reduce(ρ) + p = sortcontentsperm(Λ) + Q = Q[:,p] + Λ = Λ[p,:] n = length(ρ.generators)+1 Q̃ = similar(Q) @@ -386,20 +480,18 @@ function blockdiagonalize(ρ::Representation) c = contents2partition(Λ) k = 0 - for pⱼ in partitions(n) + for pⱼ in union(c) j = findall(isequal(pⱼ), c) - if !isempty(j) - Qⱼ = Q[:,j] - ρⱼ = Representation(map(g -> Qⱼ'*g*Qⱼ, ρ.generators)) - Q̃ⱼ = singlemultreduce(ρⱼ) - m = length(j) - Q̃[:,k+1:k+m] = Qⱼ * Q̃ⱼ - irrep = Representation(pⱼ) - for ℓ = 1:n-1 - ρd[ℓ][k+1:k+m,k+1:k+m] = blockdiag(fill(irrep.generators[ℓ], m÷size(irrep,1))...) - end - k += m + Qⱼ = Q[:,j] + ρⱼ = conjugate(ρ, Qⱼ) + Q̃ⱼ = singlemultreduce_blockdiag(ρⱼ, Representation(pⱼ)) + m = length(j) + Q̃[:,k+1:k+m] = Qⱼ * Q̃ⱼ + irrep = Representation(pⱼ) + for ℓ = 1:n-1 + ρd[ℓ][k+1:k+m,k+1:k+m] = blockdiag(fill(irrep.generators[ℓ], m÷size(irrep,1))...) end + k += m end Representation(ρd), Q̃ end diff --git a/test/runtests.jl b/test/runtests.jl index 81d7e07..202dada 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ -using NumericalRepresentationTheory, Permutations, LinearAlgebra, SparseArrays, Test -import NumericalRepresentationTheory: gelfandbasis, canonicalprojection +using NumericalRepresentationTheory, Permutations, LinearAlgebra, SparseArrays, BlockBandedMatrices, Test +import NumericalRepresentationTheory: gelfandbasis, canonicalprojection, singlemultreduce, singlemultreduce_blockdiag, conjugate, gelfand_reduce, contents2partition, sortcontentsperm @testset "Representations" begin σ = Partition([3,3,2,1]) @@ -28,7 +28,7 @@ import NumericalRepresentationTheory: gelfandbasis, canonicalprojection @testset "Rotate irrep" begin ρ = Representation(3,2,1,1) λ,Q = blockdiagonalize(ρ) - @test Q ≈ I + @test Q ≈ -I || Q ≈ I @test multiplicities(ρ)[Partition(3,2,1,1)] == 1 Q = qr(randn(size(ρ,1), size(ρ,1))).Q @@ -42,6 +42,39 @@ import NumericalRepresentationTheory: gelfandbasis, canonicalprojection g = Permutation([1,2,3]) @test ρ(g) == I(2) end + + @testset "group by rep" begin + s = standardrepresentation(4) + ρ = s ⊗ s + Λ = gelfand_reduce(ρ)[1] + p = sortcontentsperm(Λ) + @test contents2partition(Λ[p,:]) == [fill(Partition(2,1,1),3); fill(Partition(2,2),2); fill(Partition(3,1),9); fill(Partition(4),2)] + end + + @testset "singlemultreduce_blockdiag" begin + σ = Representation(3,2,1) + ρ = blockdiag(σ, σ) + Q = qr(randn(size(ρ))).Q + ρ = conjugate(ρ, Q) + Λ, Q̃ = gelfand_reduce(ρ) + p = sortcontentsperm(Λ) + ρ = conjugate(ρ, Q̃[:,p]) + Q = singlemultreduce(ρ) + + Q̃ = singlemultreduce_blockdiag(ρ, σ) + @test Q̃'Q̃ ≈ I + for (σ_k,g) in zip(blockdiag(σ,σ).generators,ρ.generators) + @test Q̃'g*Q̃ ≈ Q'g*Q ≈ σ_k + end + + ρ = blockdiag(σ, σ) + Q = qr(randn(size(ρ))).Q + ρ = conjugate(ρ, Q) + λ,Q = blockdiagonalize(ρ) + for (σ_k,g) in zip(blockdiag(σ,σ).generators,ρ.generators) + @test Q'g*Q ≈ σ_k + end + end end @testset "Canonical Projection" begin @@ -62,9 +95,10 @@ end Q = [Q_1 Q_2]' @test Q'Q ≈ I ρ_2 = Representation(3,1) - @test Q*ρ.generators[1]*Q' ≈ blockdiag(sparse(I(1)),ρ_2.generators[3]) - @test Q*ρ.generators[2]*Q' ≈ blockdiag(sparse(I(1)),ρ_2.generators[2]) - @test Q*ρ.generators[3]*Q' ≈ blockdiag(sparse(I(1)),ρ_2.generators[1]) + # we don't guarantee the ordering + @test_broken Q*ρ.generators[1]*Q' ≈ blockdiag(sparse(I(1)),ρ_2.generators[3]) + @test_broken Q*ρ.generators[2]*Q' ≈ blockdiag(sparse(I(1)),ρ_2.generators[2]) + @test_broken Q*ρ.generators[3]*Q' ≈ blockdiag(sparse(I(1)),ρ_2.generators[1]) end @testset "tensor" begin s = standardrepresentation(4) @@ -84,16 +118,16 @@ end P_211 = canonicalprojection(Partition(2,1,1), ρ) @test P_211^2 ≈ P_211 - @test P_4*Q[:,1:2] ≈ Q[:,1:2] - @test norm(P_4*Q[:,3:end]) ≤ 1E-15 - @test norm(P_31*Q[:,1:2]) ≤ 1E-15 - @test P_31*Q[:,3:3 +8] ≈ Q[:,3:3 +8] - @test norm(P_31*Q[:,3+9:end]) ≤ 1E-15 - @test norm(P_22*Q[:,1:11]) ≤ 1E-15 - @test P_22*Q[:,12:13] ≈ Q[:,12:13] - @test norm(P_22*Q[:,14:end]) ≤ 1E-14 - @test norm(P_211*Q[:,1:13]) ≤ 1E-14 - @test P_211*Q[:,14:end] ≈ Q[:,14:end] + @test P_4*Q[:,end-1:end] ≈ Q[:,end-1:end] + @test norm(P_4*Q[:,1:end-2]) ≤ 1E-15 + @test norm(P_31*Q[:,end-1:end]) ≤ 1E-15 + @test P_31*Q[:,end-10:end-2] ≈ Q[:,end-10:end-2] + @test norm(P_31*Q[:,1:end-11]) ≤ 1E-15 + @test norm(P_22*Q[:,1:3]) ≤ 1E-14 + @test P_22*Q[:,4:5] ≈ Q[:,4:5] + @test norm(P_22*Q[:,6:end]) ≤ 1E-14 + @test norm(P_211*Q[:,4:end]) ≤ 1E-14 + @test P_211*Q[:,1:3] ≈ Q[:,1:3] Q_31 = svd(P_31).U[:,1:9] @@ -118,9 +152,9 @@ end @test Q̃_31'Q̃_31 ≈ I # Q̃_31 spans same space as columns of Q - @test norm(Q̃_31'Q[:,1:2]) ≤ 1E-14 - @test rank(Q̃_31'Q[:,3:3 +8]) == 9 - @test norm(Q̃_31'Q[:,12:end]) ≤ 1E-14 + @test norm(Q̃_31'Q[:,end-1:end]) ≤ 1E-14 + @test rank(Q̃_31'Q[:,end-10:end-2]) == 9 + @test norm(Q̃_31'Q[:,1:end-11]) ≤ 1E-14 # this shows the action of ρ acts on each column space separately @test Q̃_31'ρ.generators[1]*Q̃_31 ≈ Diagonal(Q̃_31'ρ.generators[1]*Q̃_31) @@ -132,9 +166,9 @@ end Q̃[:,3:3 +8] = Q̃_31 ρ_22 = Representation(2,2) ρ_211 = Representation(2,1,1) - @test Q'*ρ.generators[1]*Q ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[1],3)..., ρ_22.generators[1], ρ_211.generators[1]) - @test Q'*ρ.generators[2]*Q ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[2],3)..., ρ_22.generators[2], ρ_211.generators[2]) - @test Q'*ρ.generators[3]*Q ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[3],3)..., ρ_22.generators[3], ρ_211.generators[3]) + @test_skip Q'*ρ.generators[1]*Q ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[1],3)..., ρ_22.generators[1], ρ_211.generators[1]) + @test_skip Q'*ρ.generators[2]*Q ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[2],3)..., ρ_22.generators[2], ρ_211.generators[2]) + @test_skip Q'*ρ.generators[3]*Q ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[3],3)..., ρ_22.generators[3], ρ_211.generators[3]) @test_skip Q̃'*ρ.generators[1]*Q̃ ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[1],3)..., ρ_22.generators[1], ρ_211.generators[1]) @test_skip Q̃'*ρ.generators[2]*Q̃ ≈ blockdiag(sparse(I(2)),fill(ρ_31.generators[2],3)..., ρ_22.generators[1], ρ_211.generators[1])