diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md
index 2a853124a..b91e0e68b 100644
--- a/docs/src/Changelog.md
+++ b/docs/src/Changelog.md
@@ -42,6 +42,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 ### Performance
 
 - Recoupling matrices are stored such that they can leverage BLAS, also for mixed complex tensor with real recoupling coefficients.
+- Constructing a `FusionTreeBlock` only takes the fusion product of the shorter side, and `fusionblocks` groups the cached fusion trees of the space instead of constructing a block for every combination of sectors (internal)
 
 ## [0.17.2](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.1...v0.17.2) - 2026-09-20
 
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index 7b97223ec..093d9daca 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -140,18 +140,29 @@ function FusionTreeBlock{I}(
     end
 
     # note: cs is unique due to ⊗ returning unique fusion results
-    cs = if N₁ == 0
-        sort!(collect(filter(isunit, ⊗(uncoupled[2]...))))
-    elseif N₂ == 0
-        sort!(collect(filter(isunit, ⊗(uncoupled[1]...))))
+    # a superset of the coupled sectors suffices, as there are no fusion trees for the others:
+    # the unit of the other side, or the fusion product of the shorter side, which avoids that of
+    # the longer side (`vec`: products of `ProductSector`s have a shape)
+    cs = if N₁ == 0 || N₂ == 0
+        # the only unit that can couple: on the left of the first and on the right of the last sector
+        us = N₁ == 0 ? uncoupled[2] : uncoupled[1]
+        leftunit(first(us)) == rightunit(last(us)) ? [leftunit(first(us))] : I[]
     else
-        sort!(collect(intersect(⊗(uncoupled[1]...), ⊗(uncoupled[2]...))))
+        sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...))))
     end
 
-    for c in cs,
-            f₂ in fusiontrees(uncoupled[2], c, isdual[2]),
-            f₁ in fusiontrees(uncoupled[1], c, isdual[1])
-        push!(trees, (f₁, f₂))
+    for c in cs
+        # the shorter side has fusion trees for every candidate, the longer side not necessarily
+        # (with several units, `c` must also have the units at the ends of the longer side)
+        us = N₁ <= N₂ ? uncoupled[2] : uncoupled[1]
+        UnitStyle(I) isa SimpleUnit || N₁ == 0 || N₂ == 0 ||
+            (leftunit(c) == leftunit(first(us)) && rightunit(c) == rightunit(last(us))) || continue
+        longer = N₁ <= N₂ ? fusiontrees(uncoupled[2], c, isdual[2]) : fusiontrees(uncoupled[1], c, isdual[1])
+        isempty(longer) && continue
+        for f₂ in fusiontrees(uncoupled[2], c, isdual[2]),
+                f₁ in fusiontrees(uncoupled[1], c, isdual[1])
+            push!(trees, (f₁, f₂))
+        end
     end
     return FusionTreeBlock(trees)
 end
diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl
index b41bd45be..dc5bd589b 100644
--- a/src/spaces/homspace.jl
+++ b/src/spaces/homspace.jl
@@ -210,13 +210,19 @@ grouped by their uncoupled charges.
 function fusionblocks(W::HomSpace)
     I = sectortype(W)
     N₁, N₂ = numout(W), numin(W)
-    isdual_src = (map(isdual, codomain(W)), map(isdual, domain(W)))
-    fblocks = Vector{FusionTreeBlock{I, N₁, N₂, fusiontreetype(I, N₁, N₂)}}()
-    for dom_uncoupled_src in sectors(domain(W)), cod_uncoupled_src in sectors(codomain(W))
-        fs_src = FusionTreeBlock{I}((cod_uncoupled_src, dom_uncoupled_src), isdual_src)
-        isempty(fs_src) || push!(fblocks, fs_src)
+    F = fusiontreetype(I, N₁, N₂)
+    # group the (cached) fusion trees of `W` by their uncoupled sectors
+    groups = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}()
+    for (f₁, f₂) in fusiontrees(W)
+        push!(get!(Vector{F}, groups, (f₁.uncoupled, f₂.uncoupled)), (f₁, f₂))
     end
-    return fblocks
+    # in the order of `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner), each of
+    # which runs fastest over its first space
+    position(V) = Dict{I, Int}(s => i for (i, s) in enumerate(sectors(V)))
+    codpos, dompos = map(position, codomain(W)), map(position, domain(W))
+    rank((u₁, u₂)) = (reverse(map(getindex, dompos, u₂)), reverse(map(getindex, codpos, u₁)))
+    us = collect(keys(groups))
+    return map(i -> FusionTreeBlock(groups[us[i]]), sortperm(map(rank, us))) # each rank computed once
 end
 
 function diagonalblockstructure(W::HomSpace)
diff --git a/test/symmetries/doubletree.jl b/test/symmetries/doubletree.jl
index 4fbe43e90..908e0497a 100644
--- a/test/symmetries/doubletree.jl
+++ b/test/symmetries/doubletree.jl
@@ -235,6 +235,21 @@ using TensorKitSectors
         end
     end
 
+    UnitStyle(I) isa SimpleUnit && @testset "Double fusion tree: blocks of a HomSpace" begin
+        V = Vect[I](Dict(s => 1 for s in (out..., incoming)))
+        W = V ⊗ V' ← V ⊗ V' # never empty: the unit couples both sides
+        trees = collect(fusiontrees(W))
+        blocks = TK.fusionblocks(W)
+        # every block holds the trees of `W` with its uncoupled sectors, in the same order
+        for block in blocks
+            @test fusiontrees(block) == filter(f -> (f[1].uncoupled, f[2].uncoupled) == block.uncoupled, trees)
+        end
+        @test sum(length, blocks) == length(trees)
+        # the blocks are ordered as `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner)
+        uncoupled = [(u₁, u₂) for u₂ in sectors(domain(W)) for u₁ in sectors(codomain(W))]
+        @test map(block -> block.uncoupled, blocks) == filter(u -> any(f -> (f[1].uncoupled, f[2].uncoupled) == u, trees), uncoupled)
+    end
+
     @testset "Double fusion tree: planar trace" begin
         if FusionStyle(I) isa UniqueFusion
             f1, f1 = src