Improve powerSet performance

[treeowl]

Obviously, there's only so much we can do, but we can do some. The most obvious optimization is to use a version of insertMin that takes a singleton tree as an argument instead of an element. But I can't help wondering if we can do better. The problem is obviously $\Omega(2^n)$, and our solution is $O(2^n\log n)$. Can we close that gap by either improving our solution or proving a tight(er) lower bound?

[treeowl]

@meooow25 , I started rambling about this in comments on your PR. Feel very welcome to think about it.

[treeowl]

We can certainly improve the space requirement for the result to $O(2^n)$, though I don't yet know whether or we can do so without making the time blow up. In particular, each set in the result with $k \ge 3$ elements can be built by bining two sets in the result of equal or nearly equal size with an element. So the trick is to get the right sets and the right element where we need them when we need them.

[treeowl]

Continuing to think out loud: imagine we build all the singletons first, then the 2-sets, then the 3-sets, etc. When we're working on the $k$-sets, we can hold on to the collections of $\lfloor k/2 \rfloor$-sets and $\lceil k/2 \rceil$-sets. The question is, how do we know which one to combine with which, and with which root element? I've drawn out the first few stages of this process, but it's hard to draw enough of it to look for patterns because the diagram grows exponentially! In case it helps anyone, here's a diagram for powerSet [1..5]. Each starting column represents a set size, and sets are lexicographically ordered from top to bottom.

[]
  [1]
    [1,2]
      [1,2,3]
        [1,2,3,4]
          [1,2,3,4,5]
        [1,2,3,5]
      [1,2,4]
        [1,2,4,5]
      [1,2,5]
    [1,3]
      [1,3,4]
        [1,3,4,5]
      [1,3,5]
    [1,4]
      [1,4,5]
    [1,5]
  [2]
    [2,3]
      [2,3,4]
        [2,3,4,5]
      [2,3,5]
    [2,4]
      [2,4,5]
    [2,5]
  [3]
    [3,4]
      [3,4,5]
    [3,5]
  [4]
    [4,5]
  [5]

[treeowl]

For the next ... while ... I will assume that we're taking the power set of [1..n]. That can be generalized by tracking the position (in the original set) of the last element of a subset. Suppose we're building from columns j and k (which will either be identical or adjacent), so that each element we build will look like p ++ [e] ++ q, where p is from j and q is from k. We walk p through j from its beginning up until we reach a point where the largest element of p is too big to be able to choose e and still have enough room to take q. For each p, we need to be able to find the matching qs, and we have to do it efficiently. I don't think we can afford to walk k from the start for each p. So let's walk k once to break it up into lists by first element, and put those lists in an array indexed by first element. Now we can quickly find where the qs start for any given p.

The hardest remaining part seems to be fitting the collections of elements of each size together into a single Set. There's that wobbling back and forth in size. We have to know when to zig and when to zag.

[treeowl]

@ekmett Do you have any ideas, either in the direction I'm pointing or another one?

[meooow25]

I don't fully understand your last comment, but I can follow the ones before it. The diagram is also interesting.

Here has been my line of thought:

The current algorithm isn't useful in getting to $O(2^n)$ because of insertMin. If that were $O(1)$, we would be good. So this works in optimal time for lists, for instance, but not sets.

To construct sets in $O(1)$, we need left and right subtrees that are balanced wrt each other and the root element.
We can have an array :: Array Int (Set a) where the index is a bitmask denoting indices of the set of elements in the original set. Then we can build the sets for a mask by finding the middle set bit, taking that as root, then indexing into the array with the lower and upper parts of the mask for left and right subtrees.

If we find the middle bit by looping over the positions, we get $O(n 2^n)$ time but $O(2^n)$ memory which is nice.
If we instead binary search for the middle bit, given we have something like $O(1)$ popCount, we are back to $O(2^n \log n)$ time but still with $O(2^n)$ memory.
I don't know of a way to find the middle bit in $O(1)$, if we could do that we would be golden.

[meooow25]

It's actually not difficult.

middleBit 1 = 0
middleBit n
  | even n    = middleBit (n `div` 2) + 1
  | otherwise = middleBit (n - 1)

We can memoize for $O(1)$.

So that's at least one way to build the power set in $O(2^n)$.
But I overlooked the fact that we need to convert from the mask order to the sorted order of sets. Trying to figure that out.

[treeowl]

I don't understand what you mean about a middle bit.

A bit more explanation of where I was heading: For each suitable left child, we choose one or more suitable roots. For each such root, we choose one or more right children. For instance, suppose we're walking the collection of 2-sets and looking at the collection of 3-sets to form the collection of 6-sets (we could actually do it either way, depending on whether we want our trees to be heavier on the left or on the right). If we've walked up to [1,3] in the collection of 2-sets, then we need to consider, in turn, each root 4, 5, .... We can't use too big a root, because we need to have at least three greater elements in the right child, so we can stop considering roots once they get too big. Suppose we are at root 5, so we have a left subtree of [1,3] and a root of 5. We now need to use, as right children, all the 3-trees whose first element is at least 6. My array concept was intended to make it cheap to jump to those 3-trees, though I'm not sure if it's theoretically necessary.

[treeowl]

Ah, the pattern is obvious now that I see it. Silly me! We always go up by one in length after a subset that does not include the maximum value. We always go down by one in length otherwise. So once we have all the subsets, getting them in order is actually easy. We can make them a set using fromDistinctAscList (or a version thereof that knows the size in advance). Is there a better way to form the set in the end? Maybe, but let's get to something that works first and think about constant factors (and GC effects) later.

[meooow25]

Thanks for the explanation, I can see how you want to construct the sets, but I don't understand how it all ties together. Perhaps some code/pseudocode can make things clearer.

I don't understand what you mean about a middle bit.

The middle set bit of 0b11011001 is bit 4, for example (0-indexed). To construct the set for this mask, we take element 4 as the root, the set we constructed for 0b1001 as the left subtree and that for 0b11000000 as the right subtree.

Here is my idea implemented. Not benchmarked or optimized.

import Data.Bits (bit, xor, (.&.), popCount, finiteBitSize, countLeadingZeros)
import qualified Data.Array.Unboxed as UA
import qualified Data.List as L
import Data.Set.Internal

-- O(2^n)
powerSet :: forall a. Set a -> Set (Set a)
powerSet xs0 = fromDistinctAscList (L.map (sets UA.!) sortedMasks)
  where
    n = size xs0
    xs = UA.listArray (0, n - 1) (toList xs0) :: UA.Array Int a

    -- TODO: Use unboxed array for better constant factor
    middleBit :: UA.Array Int Int
    middleBit = UA.listArray (1, bit n - 1) (L.map f [1 .. bit n - 1])
      where
        f 1 = 0
        f msk
          | even msk            = (middleBit UA.! (msk `quot` 2)) + 1
          | even (popCount msk) = middleBit UA.! (msk - 1)
          | otherwise           = let x = middleBit UA.! (msk - 1)
                                      -- return the next lowest set bit after x
                                      msk' = msk .&. (bit x - 1)
                                  in finiteBitSize msk' - 1 - countLeadingZeros msk'

    sets :: UA.Array Int (Set a)
    sets = UA.listArray (0, bit n - 1) (L.map f [0 .. bit n - 1])
      where
        f 0 = Tip
        f msk = bin x (sets UA.! mskl) (sets UA.! mskr)
          where
            m = middleBit UA.! msk
            x = xs UA.! m
            mskl = msk .&. (bit m - 1)
            mskr = msk `xor` bit m `xor` mskl

    -- TODO: Use an unboxed array for better constant factor
    sortedMasks :: [Int]
    sortedMasks = 0 : L.foldl' step [] [n-1, n-2 .. 0]
      where
        step msks i = bit i : L.map (bit i +) msks ++ msks

[jwaldmann]

It's fascinating to follow your discussion, but, forgive me, where's the use case?

Quite unscientifically, I checked the (approx. 36.000) *.hs that I got from unpacking $HOME/.cabal/packages/**/*.tar.gz (just what I happen to have lying around) - the only mention of powerSet is its implementation (and tests) in containers.

Another quick experiment shows

ghci> foldMap (Sum . length) $ S.powerSet $ S.fromList [1..20::Int]
Sum {getSum = 10485760}
(0.99 secs, 483,502,672 bytes)

ghci> foldMap (Sum . length) $ L.subsequences [1..20::Int]
Sum {getSum = 10485760}
(0.29 secs, 278,620,464 bytes)

I am reading this as: with current implementation, tree balancing doubles space, and triples time, w.r.t. no balancing (and no trees at all). And that's fine?

We really can't apply powerSet to anything substantially larger, because then it wouldn't fit in memory? (But it has to, since the outer tree, and all inner trees, are strict?)

Perhaps something Set a -> [Set a] (producing a lazy list) is more useful? The first candidate is map S.fromList . L.subsequences. (And their "maximally lazy" order is nice.)

Well, there's always the application as an exercise in a data structures course, so by all means, go ahead...

[meooow25]

I agree that there's not much utility to this, personally I've just been treating this as a puzzle 😁

We really can't apply powerSet to anything substantially larger, because then it wouldn't fit in memory?

Probably so in practice. The formal limit for the input set size can be considered to be word_size - 1 since the output size must not exceed maxBound :: Int.

[treeowl]

@jwaldmann This is absolutely just a puzzle. IIRC, I added powerSet to Data.Set because @ekmett thought it should be there, and offered the simple implementation we currently use. Something Set -> [Set] would certainly be accepted if there's a use for it and a good implementation. Since it's not a "fundamental operation on mathematical sets", I don't think it can slip in under the door like powerSet.

Going back to powerSet: Dropping from $O(2^n \log n)$ space to $O(2^n)$ space (with optimal sharing) should buy usability for a few extra input elements. On the time side, the logarithmic factor probably doesn't matter, so if someone has an $O(2^n)$ space/ $O(2^n \log n)$ time algorithm that actually runs faster than whatever $O(2^n)$ time algorithm we find, we can use the "theoretically slower" one instead.

[treeowl]

@meooow25 Your implementation looks interesting, but I can't make head or tail of how it works. Could you explain? How does it achieve lexicographical order? I also noticed a type error: unboxed arrays can't hold Sets.