16

[This is related to Minimum set cover ]

I would like to solve the following puzzle by computer for small size of n. Consider all 2^n binary vectors of length n. For each one you delete exactly n/3 of the bits, leaving a binary vector length 2n/3 (assume n is an integer multiple of 3). The goal is to choose the bits you delete so as to minimize the number of different binary vectors of length 2n/3 that remain at the end.

For example, for n = 3 the optimal answer is 2 different vectors 11 and 00. For n = 6 it is 4, for n = 9 it is 6 and for n = 12 it is 10.

I had previously attempted to solve this problem as a minimum set cover problem of the following sort. All the lists contain only 1s and 0s.

I say that a list A covers a list B if you can make B from A by inserting exactly x symbols.

Consider all 2^n lists of 1s and 0s of length n and set x = n/3. I would like to compute a minimal set of lists of length 2n/3 that covers them all. David Eisenstat provided code that converted this minimal set cover problem into a Mixed Integer Programming Problem that could be fed into CPLEX (or http://scip.zib.de/ which is open source).

from collections import defaultdict
from itertools import product, combinations

def all_fill(source, num):
    output_len = (len(source) + num)
    for where in combinations(range(output_len), len(source)):
        poss = ([[0, 1]] * output_len)
        for (w, s) in zip(where, source):
            poss[w] = [s]
        for tup in product(*poss):
            (yield tup)

def variable_name(seq):
    return ('x' + ''.join((str(s) for s in seq)))
n = 12
shortn = ((2 * n) // 3)
x = (n // 3)
all_seqs = list(product([0, 1], repeat=shortn))
hit_sets = defaultdict(set)
for seq in all_seqs:
    for fill in all_fill(seq, x):
        hit_sets[fill].add(seq)
print('Minimize')
print(' + '.join((variable_name(seq) for seq in all_seqs)))
print('Subject To')
for (fill, seqs) in hit_sets.items():
    print(' + '.join((variable_name(seq) for seq in seqs)), '>=', 1)
print('Binary')
for seq in all_seqs:
    print(variable_name(seq))
print('End')

The problem is that if you set n=15 then the instance it outputs is too large for any solver I can find. Is there a more efficient way of solving this problem so I can solve n=15 or even n = 18?

Community
  • 1
  • 1
Simd
  • 19,447
  • 42
  • 136
  • 271
  • 2
    I don't get it, for n=3, no matter which bit you delete, you remain with the four vectors 00 01 10 11, what am I getting wrong? – Ron Teller Oct 11 '13 at 09:29
  • @RonTeller Take all 8 vectors of length 3. For 000, delete one 0 bit. For every vector with exactly one bit set to 1 delete that bit. So all of those went to 00. For 111, delete one 1 bit. For every vector with exactly one bit set to 0, delete that bit. So all of those went to 11 and you are done. – Simd Oct 11 '13 at 09:32
  • @RonTeller The bits you delete from each vector can be different. Maybe that wasn't clear? – Simd Oct 11 '13 at 09:42
  • Exactly, now it's clear – Ron Teller Oct 11 '13 at 09:59
  • 5
    The underlying mathematical problem was [posted to Math.SE](http://math.stackexchange.com/q/497343/3111) last month, and when it proved to be difficult, [cross-posted to MathOverflow](http://mathoverflow.net/q/142857/10700). Also the problem seems to have been [posted to StackOverflow but deleted](http://stackoverflow.com/questions/18869889/deletion-puzzle-involving-binary-vectors). I'm okay with asking "a more efficient way of solving this problem," but it's good to have some links to previous efforts. – hardmath Oct 11 '13 at 10:25
  • @hardmath Yes although the only link that contains a previous attempt to write code for this problem is the one I gave at the top of the question. – Simd Oct 11 '13 at 10:30
  • 1
    Would code to directly solve the problem "more efficient[ly]", other than the metaprogramming approach above (writing Python that constructs CPLEX input), be of interest (provided the case n=15 is feasible)? – hardmath Oct 11 '13 at 11:11
  • @hardmath Yes absolutely! In fact that would be much better. – Simd Oct 11 '13 at 11:16
  • @felix, can you provide solution for n=9? – barmatat Oct 11 '13 at 23:24
  • @barmatat Sure. 000000, 001111, 011000, 100111, 110000, 111111 is one minimal solution is for n = 9. – Simd Oct 12 '13 at 07:45
  • @felix, this set of six vectors seems to contradict that you've written, the optimal answer for n=9 is 8 – elias Oct 18 '13 at 09:34
  • @elias why you do you say that? Which vector of length 9 can't be made into one of those 6 vectors by deleting 3 bits? – Simd Oct 18 '13 at 10:04
  • @felix quote from the original question: For example, for n = 3 the optimal answer is 2 different vectors 11 and 00. For n = 6 it is 4, for n = 9 it is 8 and for n = 12 it is 10. – elias Oct 18 '13 at 12:09
  • @elias Thanks that was a typo. Fixed. – Simd Oct 18 '13 at 12:32

2 Answers2

4

This doesn't solve your problem (well, not quickly enough), but you're not getting many ideas and someone else may find something useful to build on here.

It's a short pure Python 3 program, using backtracking search with some greedy ordering heuristics. It solves the N = 3, 6, and 9 instances very quickly. It finds a cover of size 10 for N=12 quickly too, but will apparently take a much longer time to exhaust the search space (I'm out of time for this, and it's still running). For N=15, the initialization time is already slow.

Bitstrings are represented by plain N-bit integers here, so consume little storage. That's to ease recoding in a faster language. It does make heavy use of sets of integers, but no other "advanced" data structures.

Hope this helps someone! But it's clear that the combinatorial explosion of possibilities as N increases ensures that nothing will be "fast enough" without digging deeper into the mathematics of the problem.

def dump(cover):
    for s in sorted(cover):
        print("    {:0{width}b}".format(s, width=I))

def new_best(cover):
    global best_cover, best_size
    assert len(cover) < best_size
    best_size = len(cover)
    best_cover = cover.copy()
    print("N =", N, "new best cover, size", best_size)
    dump(best_cover)

def initialize(N, X, I):
    from itertools import combinations
    # Map a "wide" (length N) bitstring to the set of all
    # "narrow" (length I) bitstrings that generate it.
    w2n = [set() for _ in range(2**N)]
    # Map a narrow bitstring to all the wide bitstrings
    # it generates.
    n2w = [set() for _ in range(2**I)]
    for wide, wset in enumerate(w2n):
        for t in combinations(range(N), X):
            narrow = wide
            for i in reversed(t):  # largest i to smallest
                hi, lo = divmod(narrow, 1 << i)
                narrow = ((hi >> 1) << i) | lo
            wset.add(narrow)
            n2w[narrow].add(wide)
    return w2n, n2w

def solve(needed, cover):
    if len(cover) >= best_size:
        return
    if not needed:
        new_best(cover)
        return
    # Find something needed with minimal generating set.
    _, winner = min((len(w2n[g]), g) for g in needed)
    # And order its generators by how much reduction they make
    # to `needed`.
    for g in sorted(w2n[winner],
                    key=lambda g: len(needed & n2w[g]),
                    reverse=True):
        cover.add(g)
        solve(needed - n2w[g], cover)
        cover.remove(g)

N = 9  # CHANGE THIS TO WHAT YOU WANT

assert N % 3 == 0
X = N // 3 # number of bits to exclude
I = N - X  # number of bits to include

print("initializing")
w2n, n2w = initialize(N, X, I)
best_cover = None
best_size = 2**I + 1  # "infinity"
print("solving")
solve(set(range(2**N)), set())

Example output for N=9:

initializing
solving
N = 9 new best cover, size 6
    000000
    000111
    001100
    110011
    111000
    111111

Followup

For N=12 this eventually finished, confirming that the minimal covering set contains 10 elements (which it found very soon at the start). I didn't time it, but it took at least 5 hours.

Why's that? Because it's close to brain-dead ;-) A completely naive search would try all subsets of the 256 8-bit short strings. There are 2**256 such subsets, about 1.2e77 - it wouldn't finish in the expected lifetime of the universe ;-)

The ordering gimmicks here first detect that the "all 0" and "all 1" short strings must be in any covering set, so pick them. That leaves us looking at "only" the 254 remaining short strings. Then the greedy "pick an element that covers the most" strategy very quickly finds a covering set with 11 total elements, and shortly thereafter a covering with 10 elements. That happens to be optimal, but it takes a long time to exhaust all other possibilities.

At this point, any attempt at a covering set that reaches 10 elements is aborted (it can't possibly be smaller than 10 elements then!). If that were done wholly naively too, it would need to try adding (to the "all 0" and "all 1" strings) all 8-element subsets of the 254 remaining, and 254-choose-8 is about 3.8e14. Very much smaller than 1.2e77 - but still way too large to be practical. It's an interesting exercise to understand how the code manages to do so much better than that. Hint: it has a lot to do with the data in this problem.

Industrial-strength solvers are incomparably more sophisticated and complex. I was pleasantly surprised at how well this simple little program did on the smaller problem instances! It got lucky.

But for N=15 this simple approach is hopeless. It quickly finds a cover with 18 elements, but makes no more visible progress for at least hours. Internally, it's still working with needed sets containing hundreds (even thousands) of elements, which makes the body of solve() quite expensive. It still has 2**10 - 2 = 1022 short strings to consider, and 1022-choose-16 is about 6e34. I don't expect it would visibly help even if this code were sped by a factor of a million.

It was fun to try, though :-)

And a small rewrite

This version runs at least 6 times faster on a full N=12 run, simply by cutting off futile searches one level earlier. Also speeds initialization, and cuts memory use by changing the 2**N w2n sets into lists (no set operations are used on those). It's still hopeless for N=15, though :-(

def dump(cover):
    for s in sorted(cover):
        print("    {:0{width}b}".format(s, width=I))

def new_best(cover):
    global best_cover, best_size
    assert len(cover) < best_size
    best_size = len(cover)
    best_cover = cover.copy()
    print("N =", N, "new best cover, size", best_size)
    dump(best_cover)

def initialize(N, X, I):
    from itertools import combinations
    # Map a "wide" (length N) bitstring to the set of all
    # "narrow" (length I) bitstrings that generate it.
    w2n = [set() for _ in range(2**N)]
    # Map a narrow bitstring to all the wide bitstrings
    # it generates.
    n2w = [set() for _ in range(2**I)]
    # mask[i] is a string of i 1-bits
    mask = [2**i - 1 for i in range(N)]
    for t in combinations(range(N), X):
        t = t[::-1]  # largest i to smallest
        for wide, wset in enumerate(w2n):
            narrow = wide
            for i in t:  # delete bit 2**i
                narrow = ((narrow >> (i+1)) << i) | (narrow & mask[i])
            wset.add(narrow)
            n2w[narrow].add(wide)
    # release some space
    for i, s in enumerate(w2n):
        w2n[i] = list(s)
    return w2n, n2w

def solve(needed, cover):
    if not needed:
        if len(cover) < best_size:
            new_best(cover)
        return
    if len(cover) >= best_size - 1:
        # can't possibly be extended to a cover < best_size
        return
    # Find something needed with minimal generating set.
    _, winner = min((len(w2n[g]), g) for g in needed)
    # And order its generators by how much reduction they make
    # to `needed`.
    for g in sorted(w2n[winner],
                    key=lambda g: len(needed & n2w[g]),
                    reverse=True):
        cover.add(g)
        solve(needed - n2w[g], cover)
        cover.remove(g)

N = 9  # CHANGE THIS TO WHAT YOU WANT

assert N % 3 == 0
X = N // 3 # number of bits to exclude
I = N - X  # number of bits to include

print("initializing")
w2n, n2w = initialize(N, X, I)

best_cover = None
best_size = 2**I + 1  # "infinity"
print("solving")
solve(set(range(2**N)), set())

print("best for N =", N, "has size", best_size)
dump(best_cover)
Tim Peters
  • 67,464
  • 13
  • 126
  • 132
  • 1
    Thank you. It's only the print commands that make it python3 as far as I can tell so pypy also works with it. This is a very nice start. For n = 15 it finds a cover of size 18 quickly with pypy which isn't bad but then gets stuck there. All 0s and all 1s are always in the solution that could be one (very) small speed up. – Simd Oct 18 '13 at 18:15
  • 1
    Right, no changes are really needed for Python 2.7.5 - some of the output just "looks funny" then, because of the `print` instances. The greedy ordering heuristics give it a good chance of finding a decent cover quickly, but it can still take forever ;-) to run to exhaustion. – Tim Peters Oct 18 '13 at 18:21
  • 1
    @felix, the greedy ordering heuristics already guarantee that "all 0" and "all 1" are the first two elements added to the cover - so there's really nothing to be gained there. – Tim Peters Oct 18 '13 at 18:34
  • In case it helps somehow, 0000001111,0000111111,0001101100,0011100011,0011111000,0100111010,0110001‌001,01‌11001110,1000010000,1001100111,1100011100,1101000110,1110000011,1111‌100000,11111‌​11000 plus all 0s and 1s is a size 17 solution your code does not find for n=15. – Simd Oct 21 '13 at 19:24
  • Also, your new code took only 4 minutes with pypy for n=12 to complete which I find very impressive. – Simd Oct 21 '13 at 19:35
  • @felix, yes, at 12 it's still lucky ;-) For n=15, I'm pretty sure it *will* find the 17-element cover, but maybe not in our lifetimes :-( Do you know yet whether 17 is optimal? – Tim Peters Oct 21 '13 at 20:49
  • I am afraid I have no real idea if 17 is optimal although my hunch is that it is. – Simd Oct 21 '13 at 22:18
  • It's been some years ;-) There is a size-16 cover for n=15: 0000011111,0001100001,0001110110,0011111100,0110000110,0110111000,0111100111 plus the bit complements of those plus all-0 and all-1. There is no smaller cover closed under bit complement. Found by a very similar program restricted to trying covers closed under bit complement, and using "limited discrepancy search" to greatly cut the time to make progress. – Tim Peters Oct 06 '20 at 16:19
-1

First consider if you have 6 bits. You can throw away 2 bits. Therefore, any pattern balance of 6-0, 5-1 or 4-2 can be converted to 0000 or 1111. In the case a 3-3 zero-one balance any pattern can be converted to one of four cases: 1000, 0001, 0111, or 1110. Therefore, one possible minimum set for 6 bits is:

0000
0001
0111
1110
1000
1111

Now consider 9 bits with 3 thrown away. You have the following set of 14 master patterns:

000000
100000
000001
010000
000010
110000
000011
001111
111100
101111
111101
011111
111110
111111

In other words, each pattern set has ones/zeros in the center, with every permutation of n/3-1 bits on each end. For example, if you have 24 bits then you will have 17 bits in the center and 7 bits on the ends. Since 2^7 = 128 you will have 4 x 128 - 2 = 510 possible patterns.

To find correct deletions there are various algorithms. One method is to find the edit distance between the current bit set and each master pattern. The pattern with the minimum edit distance is the one to convert to. This method uses dynamic programming. Another method would be to do a tree search through the patterns using a set of rules to find the matching pattern.

Tyler Durden
  • 11,156
  • 9
  • 64
  • 126
  • 1
    I am not sure this is right. For n = 6 a minimal size answer should have size 4 (not 6) with an example solution being 0000, 0111,1000,1111. I don't know what you mean in your example by a minimum set for 6 bits that has 6 vectors in it. – Simd Oct 11 '13 at 22:11
  • Yes, I guess thats true. – Tyler Durden Oct 12 '13 at 20:52