Faster iterator for planar set partitions

[tscrim]

Right now, we iterate through all planar set partitions in algebras.PlanarPartition by filtering out the non-planar diagrams. However, this is very inefficient for large values of n. We implement a recursive algorithm that works by simply taking the part {a, b, c, ...} that contains the largest element and uses the fact that we can form the planar partition by combining the planar set partitions on the remaining sets {1, ..., a-1}, {a+1, ..., b-1}, ..., which are all independent.

CC: @srdoty @fchapoton @zabrocki @anneschilling @saliola

Component: combinatorics

Keywords: set partition, diagram algebra

Author: Travis Scrimshaw

Branch/Commit: 22d8a70

Reviewer: Frédéric Chapoton

Issue created by migration from https://trac.sagemath.org/ticket/34579

[tscrim]

comment:1

sage: %timeit list(planar_diagrams_new(1))
2.24 µs ± 21.5 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
sage: %timeit list(planar_diagrams_new(2))
54.2 µs ± 279 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
sage: %timeit list(planar_diagrams_new(3))
590 µs ± 13.3 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit list(planar_diagrams_new(4))
6.05 ms ± 62.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
sage: %timeit list(planar_diagrams_new(5))
69.5 ms ± 1.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

versus old

sage: %timeit list(planar_diagrams(1))
6.36 µs ± 38.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
sage: %timeit list(planar_diagrams(2))
55.4 µs ± 333 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
sage: %timeit list(planar_diagrams(3))
1.12 ms ± 6.76 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
sage: %timeit list(planar_diagrams(4))
28.1 ms ± 169 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
sage: %timeit list(planar_diagrams(5))
872 ms ± 5.12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

So we already see a 10x speedup on n=5 (which means planar set partitions on 2n = 10 elements).

The downsides are that it no longer iterates through the set partitions in order of the number of parts and it is limited by the Python recursion depth.

A similar change could be done with the Temperley-Lieb diagrams by using the DyckWords iterator an to_noncrossing_set_partition() method, and I expect to lead to a speedup as well.

---

New commits:

22d8a70

Implement recursive algorithm for computing all planar partitions.

[tscrim]

Description changed:

--- 
+++ 
@@ -1 +1 @@
-Right now, we iterate through all planar set partitions in `algebras.PlanarPartition` by filtering out the non-planar diagrams. However, this is very inefficient for large values of `n`. We implement a recursive algorithm that works by simply taking the part `{a, b, c, ...}` that contains the largest element and uses the fact that we can form the planar partition by combining the planar set partitions on the remaining sets `{1, ..., a-1}, `{a+1, ..., b-1}`, `...`, which are all independent.
+Right now, we iterate through all planar set partitions in `algebras.PlanarPartition` by filtering out the non-planar diagrams. However, this is very inefficient for large values of `n`. We implement a recursive algorithm that works by simply taking the part `{a, b, c, ...}` that contains the largest element and uses the fact that we can form the planar partition by combining the planar set partitions on the remaining sets `{1, ..., a-1}`, `{a+1, ..., b-1}`, `...`, which are all independent.

[tscrim]

Branch: public/combinat/faster_planar_partition_iter-34579

[tscrim]

Commit: 22d8a70

[fchapoton]

comment:2

Looks good, but the definition of "planar set partition" is not clear to me. How is this different from noncrossing partitions ?

[tscrim]

comment:3

They are not different; just an alternative name.

[dcoudert]

comment:4

Could it be interesting to use some cache for recursive calls on the same sets ?

[fchapoton]

comment:5

If this is really the same as noncrossing partitions, it only depends on n, up to the unique increasing relabeling.

[tscrim]

comment:6

That is a good idea to increase the speed, but it will significantly increase the memory usage. Right now, I think this only needs to keep roughly the current planar set partition in memory. Your proposal will need to keep all k < n set partitions of n in memory.

[fchapoton]

comment:7

Just some not-well-cooked ideas:

• would it make sense to use nauty for this, if it can provide ?

• you are right about memory. Even if there is only one cached set for each integer n, the size of this set is growing like ~4^n.

• it seems to me that iterating through Dyck paths and using to_noncrossing_partition is not as fast as what you propose here

But a positive review in the current state in also possible, as it is an improvement, of course.

[tscrim]

comment:8

Replying to Frédéric Chapoton:

• would it make sense to use nauty for this, if it can provide ?

I am not sure what aspects you want to use. If it is enumerating planar graphs, this is done by plantri for connected graphs, which are all the same set partition. So I don't see how to get all set partitions from this.

• it seems to me that iterating through Dyck paths and using to_noncrossing_partition is not as fast as what you propose here

The to_noncrossing_partition() can be made a lot better. Even at its current form, it to be much slower than filtering out all Brauer diagrams by checking planarity for any n > 7:

sage: TL = algebras.TemperleyLieb(7, QQ.one())
sage: %time L = [d for d in TL.basis().keys()]
CPU times: user 572 ms, sys: 41 µs, total: 572 ms
Wall time: 571 ms
sage: len(L)
429
sage: %time L = [d.to_noncrossing_partition() for d in DyckWords(7)]
CPU times: user 10.7 ms, sys: 0 ns, total: 10.7 ms
Wall time: 10.6 ms
sage: DyckWords(7).cardinality()
429

sage: TL = algebras.TemperleyLieb(8, QQ.one())
sage: %time L = [d for d in TL.basis().keys()]
CPU times: user 8.33 s, sys: 0 ns, total: 8.33 s
Wall time: 8.33 s
sage: len(L)
1430
sage: %time L = [d.to_noncrossing_partition() for d in DyckWords(8)]
CPU times: user 29.3 ms, sys: 0 ns, total: 29.3 ms
Wall time: 29.1 ms
sage: DyckWords(8).cardinality()
1430

Some small tweaks would be made (separating out the core part of the algorithm into Cython to have it return tuples of tuples), which would likely make the speed disparity even greater.

[fchapoton]

comment:9

ok, let's move forward. I am setting the branch here to positive

[fchapoton]

Reviewer: Frédéric Chapoton

[tscrim]

comment:10

Thank you!

[vbraun]

Changed branch from public/combinat/faster_planar_partition_iter-34579 to 22d8a70