Construct fusion tree blocks from the shorter side, and group the cached trees in fusionblocks

[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]

Two notes here:

1. I think the SimpleUnit check is actually not required anymore, since we are now guaranteeing that the coupled sectors will belong to exactly one color.
2. My original reservations about the validity of using only one side for the coupled sectors is wrong here, because what is actually going on is that you moved the filter call into the for c in cs loop below, therefore getting equivalent results.

I think this does make sense, because you no longer need to generate all fusiontrees and instead can simply only verify the candidate cs the other side already restricted to.

[leburgel]

I removed the SimpleUnit check. For N₁ == 0 or N₂ == 0 the candidates are now allunits(I), otherwise the fusion product of the shorter side, for every UnitStyle. Only the longer side is checked for trees, since the shorter one has trees for every candidate.

[lkdvos]

In principle leftunit/rightunit of the other side is exactly the one you want, because I'm not sure if otherwise the fusiontree iterator might error

[lkdvos]

TODO: check that our iterator has a fast path for isempty, or is actually calling a fast path here.

[leburgel]

There doesn't seem to be a fast path, this just goes through iterate. So if there's no tree for a given coupled charge, it still goes through the intermediate channels. Do you think it's best to add a fast path for this?

[lkdvos]

I'm not entirely sure, I don't really have a concrete idea of what I was after, I just wanted to double-check that this is not doing something silly

[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?

[lkdvos]

Do we have an idea how much time this takes? If it is too much, there might be a point in caching the combination of fusiontrees and fusionblocks for non-abelian tensors

[leburgel]

After c643ae2: 0.43 s (fℤ₂ × U(1) × SU(3)) and 4.15 s (fℤ₂ × U(1) × SU(2) × U(1)) for all spaces of one CTMRG iteration. With #580 that is 4–8 % of building the braid transformers, and it runs before the threaded part. One iteration calls fusionblocks 550 times on 396 distinct spaces, so caching alone saves about 30 %. Caching it together with sectorstructure would also give the tree indices of each block, which saves the gettoken lookups of the source trees in GenericTreeTransformer.

Do you think it's worth looking into that for a follow-up?

[lkdvos]

It might actually be worth it, and feels a bit more uniform since very often fusionblocks is actually what you want for the non-abelian things while we are only caching fusiontrees