Given a list of 250 words of 4 letters each, what is the fastest way, in Python, to find a subsample of 100 words (drawn without replacement) so that the distribution of letters across the whole sample is most even?

List of words: http://pastebin.com/2kpiNAcC

Note there are only 18 different characters, and some occur much more frequently than others.

(Currently, "even letter distribution" is operationalized as the smallest difference between max and min letter counts and the smallest standard deviation of letter counts, although I'm open to suggestions.)

I've tried brute forcing it by drawing random samples, and after various refinements have arrived at the code below. But 1. I assume my brute forcing isn't optimal and 2. there probably is a smarter way, but I can't think of any smart approach.

I am interested in finding a better subsample, so either making the brute-forcing more efficient, or using a smarter algorithm could help. Not sure if this is better suited for Code Review (optimizing the brute forcing) or a site more suited for discussing a smarter algorithm.

import random

import numpy as np
import pandas as pd

# currently unused, as iterating over all combinations seems impossible
from itertools import combinations

from collections import Counter
from joblib import Parallel, delayed

#this just reads the file so "words" is a list of strings 
#corresponding to the 250 words
fname = "file.txt"
df = pd.DataFrame.from_csv(fname, sep="\t")
words = list(df.index)

# define function to randomly sample words and rate them
def sampl(r):
    s, mm = 1000000, 1000000  # arbitrary constants
    for w in range(r):
        w = random.sample(words, 100)
        counts = [v for k, v in Counter("".join(w)).items()]
        s_new = np.std(counts)
        mm_new = np.max(counts) - np.min(counts)
        # compare current subsample to benchmark within this run
        if s_new < s: 
            w_candidate_1 = w
            s = s_new
        if mm_new < mm: 
            w_candidate_2 = w
            mm = mm_new
        if s_new <= s and mm_new <= mm: 
            w_candidate_3 = w
    return (w_candidate_1, w_candidate_2, w_candidate_3)

# run the function a bunch (evaluating the results is fairly easy)
runs = 100
n_per_run = 999999
outs = Parallel(n_jobs=-10)(delayed(sampl)(r) 
                for r in [n_per_run] * runs)

# to find the best candidate within the output of the parallel job ...
s, mm = 1000000, 1000000
for w in outs[0]:
    counts = [v for k, v in Counter("".join(w)).items()]
    s_new = np.std(counts)
    mm_new = np.max(counts) - np.min(counts)
    if s_new < s: 
        w_candidate_1 = w
        s = s_new
    if mm_new < mm: 
        w_candidate_2 = w
        mm = mm_new
    if s_new <= s and mm_new <= mm: 
        w_candidate_3 = w

for ws in [w_candidate_1, w_candidate_2, w_candidate_3]:
    print([v for k, v in Counter("".join(ws)).items()])

Note: we have 10+ cores and 200+ GB RAM available to burn on this, so embarrassingly parallel approaches are useful.

  • Shouldn't the distribution of letter counts in the sample match that of the population you're sampling from? What is the purpose of your code? – jonrsharpe May 12 '15 at 9:44
  • @jonrsharpe: as the question says, to find a subsample "so that the distribution of letters across the whole sample is most even". For example, as in the other solutions, one that minimizes the range or variance of character counts. – jona May 12 '15 at 19:58
  • Yes, I've read the question! My point is that for most statistical uses you want the sample to have similar properties to the population it's representing. Hence my query: what is the purpose of your code? I can see what it does, and you've described it clearly, but what's it for? – jonrsharpe May 12 '15 at 20:01
  • It's not supposed to represent the, or any, population. Basically, we have started out with 250 words, but only have 100 slots available. It's experimental stimuli. – jona May 12 '15 at 20:35
up vote 12 down vote accepted

1. Analysis

In the general case, the problem of finding the sample of words with the "most even" distribution of letters is NP-hard. Here I'm considering a general instance of this problem to be:

Given an alphabet \$ Σ \$, a set \$ W \$ of 4-letter words over that alphabet, and an integer \$ n ≤ \left|W\right| \$, find the subset \$ W^* ⊆ W \$ with size \$ \left|W^*\right| = n \$ having the smallest variance of letter counts among all subsets of that size.

Here's a proof that this is NP-hard, by a reduction from EXACT COVER BY 3-SETS (that is, EXACT COVER restricted to sets of size 3, also known as "X3C").

Suppose that we have an instance of EXACT COVER BY 3-SETS in the form of a set \$ X \$ and a collection \$ S \$ of 3-element subsets of \$ X \$, and the problem is to determine if there is a subcollection \$ S^* ⊆ S \$ such that each element in \$ X \$ is contained in exactly one subset in \$ S^* \$.

This translates into the word problem as follows. Let \$ Σ = X ∪ S \$; let \$ W \$ be the set of four-letter words \$ abcd \$ such that \$ \{ a, b, c \} = d ∈ S \$, and let \$ n = { \left|X\right| \over 3 } \$. Then solve the word problem to find \$ W^* \$ with the smallest variance of letter counts. There is an exact cover if and only if the variance of the letter counts in \$ W^* \$ is zero. (Because there can only be one instance of any letter from \$ S \$, and so if the variance is zero, there must be exactly one instance of each letter from \$ X \$, and since there are \$ \left|X\right| \over 3 \$ words, each element of \$ X \$ is covered exactly once.)

So there's no efficient (polynomial time) algorithm that solves your problem. Exhaustive search is out, because \$ {250 \choose 100} \approx 6×10^{71} \$. So the best you can do is to use heuristics to find good approximate solutions. I'll describe one possibility in §4 below.

2. Review

  1. Having code at the top-level of a module makes it hard to test from the interactive interpreter, because whenever you reload the module, the code runs. It's best to guard top-level code with if __name__ == '__main__':.

  2. The comment for sampl doesn't explain the most important points about the behaviour of the function, namely: what arguments does it take? and what does it return?

  3. The number 100 is arbitrary, so it ought to be a parameter.

  4. When finding the minimum or maximum, it's best to avoid a starting point like this, but instead to rearrange the code into a generator that can then be passed to min. That way you can avoid arbitrary starting points like 1000000.

  5. Writing for w in range(r): is misleading, because w is not actually used (it's immediately overwritten by the random sample). It's conventional to write _ for unused loop variables.

  6. The code uses NumPy to do the statistics on the sample. But there are at most 400 letters in the sample, so this doesn't really gain much: NumPy only shows a benefit for large volumes of data. I think that using the standard Python functions statistics.pstdev, max and min would be adequate for this amount of data.

  7. Since the standard deviation is only used for comparison, you'd get the same results if you used the variance instead, and this would have the advantage of avoiding the square root.

  8. Instead of Counter("".join(w)), use itertools.chain.from_iterable and write Counter(chain.from_iterable(w)). This avoids concatenating the words. [Improved by Veedrac in comments.]

  9. Instead of:

    counts = [v for k, v in Counter("".join(w)).items()]
    

    write:

    counts = list(Counter(chain.from_iterable(w)).values())
    

    since you don't use the keys.

  10. I dislike the logic for selecting w_candidate_3. This is only updated if both scoring functions show an improvement for this sample, and so the "best" sample according to this logic will depend on the order in which samples are examined. This is a bad property for a measure to have. It would be better to find some numeric way to combine the two scores (for example, their sum or product), so that all samples are comparable according to this measure.

  11. The logic in sampl is complicated because you aren't sure of the right way to score the samples. In this kind of situation a better approach would be for sampl to take a scoring function as an argument. Then you could easily experiment with different kinds of score.

3. Revised code

from collections import Counter
from itertools import chain
import random
from statistics import pvariance

def letter_counts(words):
    """Generate the letter counts of the words."""
    return Counter(chain.from_iterable(words)).values()

def score_variance(words):
    """Score words according to variance of letter counts."""
    return pvariance(letter_counts(words))

def score_range(words):
    """Score words according to range of letter counts."""
    counts = list(letter_counts(words))
    return max(counts) - min(counts)

def best_sample(words, score, n=100, r=100):
    """Generate r (default 100) random samples of n (default 100) elements
    from words, and return the sample with the smallest score.

    """
    return min((random.sample(words, n) for _ in range(r)), key=score)

4. Hill climbing

A useful search technique in this kind of global optimization problem is hill climbing. The idea is to pick a random sample to start with, explore its neighbourhood by examining "nearby" samples, pick the best of the neighbours, and continue until no further improvement can be made. For example:

def best_sample_2(words, score, n=100, r=100):
    """Return a sample of n (default 100) elements from words, found by
    starting with a random sample and repeatedly hill climbing by
    generating r (default 100) neighbours and picking the one with
    smallest score, until no more progress is made.

    """
    words = set(words)
    current_sample = set(random.sample(words, n))
    current_score = score(current_sample)

    # Generate neighbouring samples (one word swapped).
    def neighbours():
        nonsample = words - current_sample
        for _ in range(r):
            s = set(current_sample)
            s.difference_update(random.sample(current_sample, 1))
            s.update(random.sample(nonsample, 1))
            yield s

    while True:
        new_sample = min(neighbours(), key=score)
        new_score = score(new_sample)
        if current_score <= new_score:
            # No neighbour is any better, so stop here.
            return current_sample
        current_sample, current_score = new_sample, new_score

Here's an example of how effective hill climbing can be in this problem. First, we'll use the original algorithm to pick the best of 10,000 samples.

>>> score_variance(best_sample(WORDS, score_variance, r=10000))
114.18300653594771

Now, using hill climbing:

>>> score_variance(best_sample_2(WORDS, score_variance))
23.241830065359476

The hill climbing approach reaches a local optimum. So usually it's a good idea to run it several times with different random starting points. For example, if I run it ten times and pick the best:

>>> f = lambda:best_sample_2(WORDS, score_variance)
>>> score_variance(min((f() for _ in range(10)), key=score_variance))
19.712418300653592

For further improvement, take a look at simulated annealing.

  • Massive improvement over my original approach, and clearing up my Python too! – jona May 12 '15 at 12:30

To emphasize a point I made to Gareth Rees, I shall present a totally-impractical polynomial time algorithm given a fixed alphabet \$A\$.

Let \$X \subseteq A^4\$ be the input set and \$T\$ be the number of elements to choose from \$X\$.

Define a bijective function \$f : \mathbb{Z}_{|A|} \rightarrow A\$.

Produce a \$(|A| + 1)\$-dimensional boolean n-cube \$\mathbf{C}\$ of side lengths \$|X| + 1\$, each element initialized to \$\text{false}\$.

Set the element in \$\mathbf{C}\$ at index \$(0, 0, ..., 0)\$ to \$\text{true}\$.

For each element \$c = (c_1, c_2, c_3, c_4) \in X\$, for each element \$c = \text{true} \in \mathbf{C}\$ element at index \$(i_1, i_2, ..., i_{|X|}, k)\$, set the element at index \$(i_1', i_2', ..., i_{|X|}', k+1)\$ to \$\text{true}\$, where

$$ i_n' = \begin{cases} i_n + 1 & f(n) \in c \\ i_n & \text{otherwise} \end{cases} $$

Find each set element in \$\textbf{C}\$ at index \$(i_1, i_2, ..., i_{|X|}, T)\$, and reduce with some constant-time reduction (in terms of \$|A|\$) to the index with minimum unevenness. If no such elements exist, halt.

This is time \$\mathcal{O}(|X|^{|A| + 2})\$ and space \$\mathcal{O}(|X|^{|A| + 1})\$.

This is evidently polynomial for fixed \$|A|\$, even if completely useless.


Gareth Rees did a great job at simplifying your original code, so I can only try to offer alternate strategies.

Here's a greedy algorithm. The strategy is to just always choose the next word that minimizes the score:

def greedy(words, score, n=100):
    words = set(words)
    chosen = set()

    def rate_word(word):
        return score(chosen | {word})

    for _ in range(n):
        best_word = min(words, key=rate_word)
        chosen.add(best_word)
        words.remove(best_word)

    return chosen

score_variance(greedy(WORDS, score_variance))
#>>> 16.950617283950617

This is completely inefficient; running score_variance each time results in \$\mathcal{O}(n^3)\$ running time. However, it's competitive with the other solutions and gives consistently better results. It's also easy to optimize if score stops being a black box.

Using this as the seed to hill climbing actually helps a little:

def hill_greedy(words, score, n=100, r=100):
    ...
    current_sample = greedy(words, score, n)
    ...

score_variance(hill_greedy(WORDS, score_variance))
#>>> 16.39506172839506

This is just Gareth Rees's code with greedy swapped in.

There are far better algorithms, but those will have to wait for some other time.

  • I had implemented a greedy algorithm very similar to yours... I used a list instead of a set, I guess you chose the set for performance improvements? Well done anyway – Caridorc May 13 '15 at 15:21
  • @Caridorc Set vs list would matter if I didn't call score \$\mathcal{O}(n^2)\$ times (and was a little more careful with it) - in this case I just used it for convenience. – Veedrac May 13 '15 at 16:27

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.