[leburgel]

This is a draft illustrating some potential performance improvements; the actual implementations are
very much up for discussion.

Same blocks, same order, less work:

• FusionTreeBlock{I}(uncoupled, isdual): for SimpleUnit sectors, the coupled sectors come from
  the fusion product of the shorter side only, instead of intersecting the products of both sides.
• fusionblocks(W): groups the cached fusiontrees(W) by uncoupled sectors, instead of constructing
  a (mostly empty) block for every combination of sectors.

New test: the blocks of a HomSpace hold exactly its trees, in order.

Cold construction of the recoupling data for one permutation of a tensor from a PEPS CTMRG step
(scripts below), main → this PR:

	main	this PR
fℤ₂ × U(1) × SU(3), 1 thread	5.29 s	4.96 s (1.07×)
fℤ₂ × U(1) × SU(3), 8 threads	1.58 s	1.42 s (1.11×)
fℤ₂ × U(1) × SU(2) × U(1), 1 thread	4.99 s	3.20 s (1.6×)
fℤ₂ × U(1) × SU(2) × U(1), 8 threads	3.45 s	2.40 s (1.4×)

SU(2) gains little from threads: Fsymbol for SU2Irrep goes through WignerSymbols, which takes
a lock on every call. Removing that lock (caching the F-symbols in TensorKitSectors) is scheduled
as a follow-up.

TensorKit main (ae8abfc), SUNRepresentations c02a3f3, TensorKitSectors v0.3.10, Julia 1.12.5; shared 80-core
machine, two interleaved rounds agree within 15 %.

fℤ₂ × U(1) × SU(2) × U(1) reproducer (julia -t8)

using TensorKit

# the spaces of one index permutation in a CTMRG step for a fℤ₂ × U(1) × SU(2) × U(1)-symmetric PEPS
const I = FermionParity ⊠ U1Irrep ⊠ SU2Irrep ⊠ U1Irrep
s(f, q, j, q′) = FermionParity(f) ⊠ U1Irrep(q) ⊠ SU2Irrep(j) ⊠ U1Irrep(q′)
V1 = Vect[I](s(0, 0, 1/2, 1) => 1, s(1, -1, 0, 0) => 1, s(0, 2, 0, 0) => 1, s(0, 0, 0, -2) => 1, s(1, 1, 1/2, -1) => 1, s(1, 1, 0, 2) => 1)
V2 = Vect[I](s(0, 0, 1/2, 1) => 1, s(0, 2, 0, 0) => 1, s(0, 0, 0, -2) => 1, s(1, 1, 1/2, -1) => 1, s(1, 1, 0, 2) => 1)
V3 = Vect[I](s(1, 1, 0, 0) => 1, s(0, 0, 1/2, -1) => 1, s(0, 0, 0, 2) => 1, s(0, 2, 1/2, 1) => 1, s(0, 2, 0, -2) => 1)
V4 = Vect[I](s(1, 1, 0, 0) => 1, s(0, 0, 1/2, -1) => 1, s(0, 0, 0, 2) => 1, s(1, -1, 1/2, 1) => 1, s(1, -1, 0, -2) => 1)
V5 = Vect[I](s(0, 0, 0, 0) => 5, s(0, 0, 1, 0) => 1, s(1, 1, 1/2, 1) => 4, s(1, -1, 1/2, -1) => 4, s(0, 0, 1/2, 3) => 1, s(1, 1, 0, -2) => 3, s(1, -1, 0, 2) => 3, s(0, 0, 1/2, -3) => 1, s(1, 1, 1, -2) => 1, s(1, -1, 1, 2) => 1, s(1, 1, 0, 4) => 1, s(1, -1, 0, -4) => 1)
V6 = Vect[I](s(0, 0, 0, 0) => 5, s(0, 0, 1, 0) => 2, s(1, 1, 1/2, 1) => 4, s(1, -1, 1/2, -1) => 3, s(0, 0, 1/2, 3) => 2, s(0, 2, 1/2, -1) => 1, s(0, -2, 1/2, 1) => 1, s(1, 1, 0, -2) => 3, s(1, -1, 0, 2) => 3, s(0, 0, 1/2, -3) => 2, s(0, 2, 0, 2) => 1, s(1, 1, 1, -2) => 1, s(0, -2, 0, -2) => 1, s(1, 1, 0, 4) => 1)
V7 = Vect[I](s(0, 0, 1/2, 1) => 1, s(1, -1, 0, 0) => 1, s(0, 0, 0, -2) => 1, s(1, 1, 1/2, -1) => 1, s(1, 1, 0, 2) => 1)
W = (V1 ⊗ V2' ⊗ V3) ← (V4' ⊗ V5 ⊗ V6' ⊗ V7)
p = ((1, 7, 4), (2, 3, 5, 6))

t = randn(W)

# cold construction of the cached recoupling data of the permutation (TensorKit's caches emptied,
# the symbols warm), with 1 and with all threads
function cold_permute(t, p, n)
    TensorKit.set_num_manipulation_threads(n)
    TensorKit.set_num_transformer_threads(n)
    TensorKit.empty_globalcaches!()
    return @elapsed permute(t, p)
end
cold_permute(t, p, Threads.nthreads()) # compilation, and the symbols
for n in unique((1, Threads.nthreads()))
    ts = [cold_permute(t, p, n) for _ in 1:3]
    println("cold permute, $n thread(s): ", round(minimum(ts); digits = 2), " s (min of 3)")
end

fℤ₂ × U(1) × SU(3) reproducer (julia -t8)

using SUNRepresentations, TensorKit
SUNRepresentations.use_disk_cache(false; persist = false)

# the spaces of one index permutation in a CTMRG step for a fℤ₂ × U(1) × SU(3)-symmetric PEPS
const I = FermionParity ⊠ U1Irrep ⊠ SU3Irrep
s(f, q, a, b) = FermionParity(f) ⊠ U1Irrep(q) ⊠ SU3Irrep(a, b)
V1 = Vect[I](s(0, 0, 1, 0) => 1, s(1, 1, 0, 1) => 1, s(1, -1, 0, 0) => 1, s(0, 2, 0, 0) => 1)
V2 = Vect[I](s(0, 0, 1, 0) => 1, s(1, -1, 0, 0) => 1, s(0, -2, 0, 1) => 1)
V3 = Vect[I](s(0, 0, 1, 0) => 1, s(1, 1, 0, 1) => 1, s(1, -1, 0, 0) => 1)
V4 = Vect[I](s(0, 0, 0, 0) => 1, s(1, -1, 0, 1) => 1, s(0, -2, 1, 0) => 1)
V5 = Vect[I](s(0, 0, 0, 0) => 2, s(0, 0, 1, 1) => 2, s(1, 1, 1, 0) => 2, s(1, -1, 0, 1) => 2, s(1, 1, 0, 2) => 1, s(0, 2, 0, 1) => 1, s(0, -2, 1, 0) => 1, s(1, -1, 2, 0) => 1)
V6 = Vect[I](s(0, 0, 0, 0) => 3, s(0, 0, 1, 1) => 2, s(1, 1, 1, 0) => 2, s(1, -1, 0, 1) => 2, s(0, 2, 0, 1) => 1, s(0, -2, 1, 0) => 1, s(0, -2, 0, 2) => 1, s(0, 2, 2, 0) => 1)
V7 = Vect[I](s(0, 0, 0, 0) => 1, s(1, 1, 1, 0) => 1, s(1, -1, 0, 1) => 1)
W = (V1 ⊗ V2 ⊗ V3) ← (V4 ⊗ V5 ⊗ V6' ⊗ V7)
p = ((1, 7, 4), (2, 3, 5, 6))

t = randn(W)

# cold construction of the cached recoupling data of the permutation (TensorKit's caches emptied,
# the F- and R-symbols warm), with 1 and with all threads
function cold_permute(t, p, n)
    TensorKit.set_num_manipulation_threads(n)
    TensorKit.set_num_transformer_threads(n)
    TensorKit.empty_globalcaches!()
    return @elapsed permute(t, p)
end
cold_permute(t, p, Threads.nthreads()) # compilation, and the symbols
for n in unique((1, Threads.nthreads()))
    ts = [cold_permute(t, p, n) for _ in 1:3]
    println("cold permute, $n thread(s): ", round(minimum(ts); digits = 2), " s (min of 3)")
end

🤖 Generated with Claude Code

[lkdvos]

Do we need to map every sector to an Int, or can we get away with just calling < on the sectors themselves since they should appear in the same order, thereby avoiding having to allocate the Dicts here?

[leburgel]

We could do this, but I checked and the isless approach is actually slower (almost a factor of 2) for the product sectors I was testing with, since that callsfindindex for each factor at every comparison. So I would be inclined to keep the Dict positions since they're faster.

Also note I fixed an issue with the sorting below, where previously the Dict lookups were run at every comparison, whereas now we run the lookup on each block once and then sort on the obtained results.

[lkdvos]

let me investigate if we can actually speed up the comparisons a bit, I think there are quite a few fast paths that might be able to get this closer together. It's also just about this looking a bit "ugly" or complicated, if we could just use comparison it's a bit easier to read

[lkdvos]

Is this specifically for the SU(3) case by the way? I would guess the findindex on the other factors is cheap (unless we are missing some fast definitions here)

[lkdvos]

QuantumKitHub/TensorKitSectors.jl#115

Mind giving it a go with that included?

[leburgel]

With QuantumKitHub/TensorKitSectors.jl#115, findindex for SU(3) goes from about 10 to 7 ns. For all spaces of one CTMRG iteration, the version with plain isless (keys computed once) takes 0.57–0.62 s instead of 0.64–0.65 s for fℤ₂ × U(1) × SU(3), and 6.3–6.9 s either way for fℤ₂ × U(1) × SU(2) × U(1); the Dict version takes 0.37–0.41 / 4.6–4.8 s with or without #115. So the remaining gap is the O(n log n) comparisons of sector tuples rather than findindex. Most of fusionblocks is grouping the trees anyway (0.27–0.33 / 3.3–3.4 s), so I'm fine with switching to the comparison version if you prefer how it reads.

[lkdvos]

I think I'd say to avoid the optimization with the integers, also because it has some implicit dependence on the spaces iterating sectors in the correct order, so if this turns out to not be dominant we should probably not. Additionally it might be worth it to check with different sector types, where it might hurt performance

[leburgel]

Switched to plain comparisons. Timings of fusionblocks over all spaces (same blocks and order in all cases):

• SU(2), A₄, Ising, Fibonacci, SU(2) ⊠ SU(2), ℤ₂ ⊠ SU(2): about the same as with the Dicts.
• SU(3) and fℤ₂ × U(1) × SU(3): about 1.8× slower (0.46 → 0.84 s for the latter, one CTMRG iteration).
• fℤ₂ × U(1) × SU(2) × U(1): also about 1.8× slower (5.6 → 10.4 s, one CTMRG iteration).

So it only shows for larger product sectors, where comparing the sector tuples is the cost. Caching the blocks (#587) recovers part of that.

[leburgel]

Small follow-up: the key of each block is now computed once (sortperm) rather than at every comparison, which takes the two product cases above from 0.84 to 0.69 s and from 10.4 to 8.1 s.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

Files with missing lines	Coverage Δ
src/fusiontrees/fusiontrees.jl	94.05% <100.00%> (-2.06%)	⬇️
src/spaces/homspace.jl	93.54% <100.00%> (+0.08%)	⬆️

... and 14 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.