Speed up getindex and findindex for ProductSector SectorValues

[borisdevos]

TLDR: it's faster.

The current Manhattan-distance enumeration was implemented just recursively, which scaled like distance^(N-1) for N sector types in the product sector. When it comes to finite product sectors within TensorKit, getindex and findindex are used frequently in constructing graded spaces.

Here, I propose a speedup that exists of two parts.

1. For a fixed Manhattan distance, there are cases where one can just combinatorically calculate the number of grid points at this distance. There's thus a shortcut now in num_manhattan_points in these cases. Otherwise, still calculate recursively. I also limited this behavior to a certain amount of product sector types just based on complexity, but it might be an unnecessary complexity (other definition) of the code.
2. Calculate the entire multi-index to Manhattan index and vice-versa just once, as it's a bijection, in such a way that still respects isless of product sectors, and cache this in a table. Then findindex and getindex are just array lookups instead of recursive calculations.

There's also a fix related to checking bounds. This is because, as of now, when getting an index outside of the length of sector values, it would get stuck in a while loop. I added checksectorindex to check the bounds. Everything's done in a way that functions can be inlined, but I don't really know how that works. I think it's worth it in this case, but please be critical about it. This infinite loop is now unreachable through the API.

A comment before the benchmarks: As of now I set the cutoff to tabulating to some value large enough that finite product sector types will most likely tabulate. Whether this tabulation matters in practice depends on Vect[I] using NTuple or SectorDict storage. Currently, every finite sectors falls within the former, but I'm looking into introducing a cutoff even in the finite case, which I hope to share soon. This is based on performance within TensorKit's functions which frequently access this storage.

Benchmarks: I checked decoding/getindex and encoding/findindex, as well as 4 functions within TensorKit which access the sector values in some way: dim, fuse, sectors and blocksectors on product spaces (not cached!). I ran this for 2 finite product sector types under the tabulation limit, an infinite product sector type, and a finite one above the limit (just the decoding/encoding).

Benchmark results

Current main

Case	decode	encode	dim(V)	sectors(V)	fuse(V,V)	blocksectors
A: Z4^⊠3(n=64)	15.52 µs	11.89 µs	5.73 µs	436 ns	2.08 ms	316 µs
B: Z8^⊠3 (n=512)	529 µs	398 µs	274 µs	2.71 µs	602 ms	169 ms
C: U1⊠SU2⊠FermionParity	61.6 µs	73.5 ns	672 ns	325 ns	372 µs	48.1 µs
D: Z8^⊠6 (n=262144)	1.76 ms	1.61 ms	/	/	/	/

Changes

Case	decode	encode	dim(V)	sectors(V)	fuse(V,V)	blocksectors
A: Z4^⊠3 (n=64)	696 ns	1.07 µs	3.91 µs	324 ns	1.71 ms	269 µs
B: Z8^⊠3 (n=512)	5.5 µs	9.04 µs	192 µs	1.66 µs	404 ms	155 ms
C: U1⊠SU2⊠FermionParity	50.0 µs	100.0 ns	588 ns	324 ns	225 µs	35.9 µs
D: Z8^⊠6 (n=262144)	15.5 ns	21.2 ns	/	/	/	/

Speedups

Case	decode	encode	dim(V)	sectors(V)	fuse(V,V)	blocksectors
A: Z4^⊠3 (n=64)	22.3×	11.2×	1.46×	1.34×	1.22×	1.18×
B: Z8^⊠3 (n=512)	96.2×	44.1×	1.42×	1.64×	1.49×	1.09×
C: U1⊠SU2⊠FermionParity	1.23×	0.74×	1.14×	1.00×	1.65×	1.34×
D: Z8^⊠6 (n=262144)	~113,000×	~76,000×	/	/	/	/

Standard errors (old / new)

Case	decode	encode	dim(V)	sectors(V)	fuse(V,V)	blocksectors
A: Z4^⊠3 (n=64)	0.6% / 0.8%	0.5% / 0.9%	0.8% / 1.0%	24.6% / 40.6%	0.3% / 2.9%	0.9% / 0.4%
B: Z8^⊠3 (n=512)	0.6% / 0.9%	0.9% / 1.1%	0.9% / 0.6%	42.9% / 0.9%	1.2% / 1.4%	0.7% / 2.5%
C: U1⊠SU2⊠FermionParity	0.7% / 1.0%	0.8% / 0.8%	1.5% / 0.8%	1.5% / 2.3%	0.5% / 0.6%	0.5% / 0.9%
D: Z8^⊠6 (n=262144)	0.8% / 1.8%	0.5% / 1.2%	/	/	/	/

Case D: Z8^⊠6 (n=262144), after execution

metric	old	new	ratio
decode, deepest index	1.76 ms	15.5 ns	~113,000×
encode, deepest index	1.61 ms	21.2 ns	~76,000×

Important observations:

• Speedup is more noticeable the larger the finite product sector type.
• For case D actually, I first had a different approach to building the table built on the recursive decoder per index. However, this was ridiculously expensive (reasoned in hindsight, of course), and added an extra factor "number of sector types in the product". build_table now enumerates the grid directly instead of this recursive decoder, which is O(n logn) due to sorting.
• The regression in case C actually doesn't matter within TensorKit, since this Manhattan code only matters when Vect[I] has NTuple storage. So this regression is only within TensorKitSectors. The last 4 columns here should really be ignored, as they don't pass the Manhattan code anyway.

I want to stress again this tabulation cutoff being arbitrary, and it really just mattering within TensorKitSectors. My preliminary benchmarks on NTuple vs SectorDict storage hint at the former being efficient from the order of ~16-32 sectors, which is way below 2^20. So really, if I end up drawing as conclusion that the cutoff should be around 16-32, then this PR is practically useful for people studying product sector types ≤ 16 sectors 🤣 Thanks for listening to my TED talk.

[codecov]

Codecov Report

❌ Patch coverage is 96.00000% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/auxiliary.jl	93.33%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/TensorKitSectors.jl	16.66% <ø> (ø)
src/product.jl	94.69% <100.00%> (+0.24%)	⬆️
src/sectors.jl	94.77% <100.00%> (+0.38%)	⬆️
src/auxiliary.jl	92.78% <93.33%> (+0.02%)	⬆️

🚀 New features to boost your workflow:

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

[lkdvos]

This is definitely cool! Do you have an idea which part is affecting the speed the most, the tabulation or the short circuit?

I also don't really know how much this actually affects performance of downstream algorithms, it seems like the affect on the functions that end up being called isn't enormous, and these seem to be very small numbers, but I don't actually know how often they are called.

For dim for example, it seems like there might be a much larger speedup from just changing the implementation here:

https://github.com/QuantumKitHub/TensorKit.jl/blob/64430c522c7da3c13c8376fc9734a8ecc054a324/src/spaces/gradedspace.jl#L92-L95

to something like:

function dim(V::GradedSpace{I, <:AbstractDict})
    init = 0 * dim(first(allunits(sectortype(V))))
    return sum((c, d) -> dim(c) * d, pairs(V.dims); init)
end
function dim(V::GradedSpace{I, <:Tuple})
    init = 0 * dim(first(allunits(sectortype(V))))
    return sum(((c, d),) -> dim(c) * d, zip(values(V). V.dims)); init)
end

In other words, completely bypassing the hash lookup that was taking place in the dim(V, c) call.

---

The main reason I'm saying this is that I'm a bit scared of adding more caches, since people tend to complain about this eating up space and this being nontransparent 🙃

[lkdvos]

Did you hit this case somewhere specifically? I did not look at this in detail, but I somehow would have expected that if the binomial function overflows also the recursive function would?

[borisdevos]

Yeah you're right, which really means that previously you could've potentially summed up incorrectly in the recursion, and now I'm catching it in a subset of cases. Although I think realistically this can't be reached, as this requires a crazy amount of product sectors, so maybe I can just remove this?

Btw I didn't hit this myself, it's something I read in the docstring of factorial so I thought I should catch it, but clearly I didn't think enough about this being useful or not.

[lkdvos]

I wonder if it might be worth it to try something similar to TensorKit's caches for this, using an LRU which is both threadsafe as well as avoids keeping too many entries around.

[lkdvos]

Do you need both here? Isn't the manhattan index 1:n, so as long as you store the multi index in sorted order by their linear index you could fold both into a single array, using the array index as the key.

[borisdevos]

I guess the tradeoffs occurring here as that you save a bit of memory removing lin (not much compared to multi, though), but looking for a multi-index given the Manhattan index requires a small binary search (which google tells me is log(n)). Is that worth it?

[lkdvos]

I also realized that I wasn't really following what was going on here anyways, and I think my suggestion is wrong and that would defeat the entire point, the ordering is important 🙃

[borisdevos]

This is definitely cool! Do you have an idea which part is affecting the speed the most, the tabulation or the short circuit?

It's the tabulation, once a sector type is tabulated you never call on the auxiliary functions again.

I also don't really know how much this actually affects performance of downstream algorithms, it seems like the affect on the functions that end up being called isn't enormous, and these seem to be very small numbers, but I don't actually know how often they are called.

For dim for example, it seems like there might be a much larger speedup from just changing the implementation here:

https://github.com/QuantumKitHub/TensorKit.jl/blob/64430c522c7da3c13c8376fc9734a8ecc054a324/src/spaces/gradedspace.jl#L92-L95

to something like: ...

In other words, completely bypassing the hash lookup that was taking place in the dim(V, c) call.

That's true, and is something I was also looking into; whether there were multiple commonly called points in TensorKit which treated ntuple and dictionary storage equally and perhaps inefficiently. However, the speedup here is already gotten purely through this tabulation. Your suggestion to dim specifically just limits the number of calls to this tabulation, so it's complementary as far as I can tell.

The main reason I'm saying this is that I'm a bit scared of adding more caches, since people tend to complain about this eating up space and this being nontransparent 🙃

Fair argument, and I think the suggestion you made concerning LRU caching is a valid approach.

[borisdevos]

Closing this as the effects of QuantumKitHub/TensorKit.jl#511, particularly the new tuple-to-sectordict cutoff, make the proposed ideas here (among others I was playing around with) completely redundant.