In our recent work, we needed to solve a maximum set coverage problem. And I found the solution so interesting that in this post, I want to discuss a very interesting algorithm to approximate solutions to a maximum set coverage problem that provide some interesting theoretical bounds. I will not go into the theoretical details that lead to the said bounds, as those are beyond the scope of this post, but you may refer to [1] and [2]. So, without further ado, let's get into it.

1   The problem of maximum coverage

The problem of maximum coverage (or max cover ) over sets is a very fundamental problem in computer science and maths with varied applications. Whether you want to compute the optimal position to place traffic sensors such that they cover the maximum possible part of the city; to find a policy that insures you against maximum financial risks; to ensure that all your code has been tested for all possible outcomes; or to find some representative objects that represent most of the dataset, you will need to solve a variant of the max cover problem.

In its purest form, the max cover problem can be defined without dealing with the earthly attachments to these applications. Therefore, we will also abandon any reference to these applications and only use some well known set theory concepts to develop our problem, and its solution.

We start with a large universe set \(U\). When we have a universe of set, we may have a family of sets \(\mathcal{P}(S)=\{S\mid S\subseteq U\}\) (note that we can only enumerate this set when \(U\) is finite, but in most applications we deal with finite sets, thus this is fine). Using the power set, we may define more families of sets such that all the set in it are subsets of the universe set, \(X\subseteq \mathcal{P}(S)\). Then we define the cover of a family of sets \(X\) called \(c(X)\) as follows.

Definition 1 (Cover of a family of sets). Let \(U\) be a set, \(X\subseteq \mathcal{P}(U)\), then \[c(X)=\left\lvert\bigcup_{S\in X} S\right\rvert\]

We say that \(X\) covers \(c(X)\) of the universe set \(U\). Visually, this may look like the schematic diagram below with the three circles representing the sets in the family \(X\).
This is a very simple concept, its just the size of unions of sets. We all learned that in high school, and have now given it a name. Nevertheless, this becomes interesting very soon because we can now define the coverage maximization problem.

Problem (Coverage maximization). Given an integer \(k\), and a universe set \(U\), and a collection of allowed sets \(A\subseteq\mathcal{P}(U)\), find a collection of sets \(S'\subseteq A\) such that \(|S'|\le k\), and the number of elements covered by \(S'\), \(c(S')\) is maximized.

This problem looks very easy, we just have to choose a set of sets such that they have the largest possible union. But, as simple as it looks, this problem is NP-hard which means that there is no polynomial time algorithm which provides a solution to this problem efficiently.

2   Greedy approximations

Even though, there do not exist any polynomial time solutions that provide the optimum results, we may still use a greedy approach to find a close approximation. How close is the approximation? In [2], the bound on the value of the greedy approximation is given as follows (under certain conditions, which do hold for max cover). \[\frac{\text{value of greedy approximation}}{\text{value of optimal solution}}\ge 1-\left(\frac{k-1}{k}\right)^k\ge \frac{e-1}{e}\]

To define the greedy algorithm, we will need to talk about the marginal gain a lot. Therefore, I define it here.

Definition 2 (Marginal cover gain). Given a family of sets \(S\), and a set \(x\subseteq \mathcal{P}(U)\), \(x\notin S\), the marginal cover gained from adding \(x\) to \(S\) is \[\sigma(S,x)=c(S\cup\{x\})-c(S).\]

The greedy algorithm amounts to iterating over the sets, computing \(\sigma(S,x)\) for each set \(x\) in each iteration, and adding the set \(x\) to \(S\) that has the maximum value of \(\sigma(S,x)\). There are some nuances, such as \(\sigma(S,x)\) should be non-zero positive, and we need to stay within the constraing of \(|S|\le k\), etc. Therefore, I define the algorithm below.

While this works great, there is still an amazing algorithmic speed-up hiding here that can be used to reduce a lot of computations and efficiently compute the greedy approximation of the maximum coverage problem.

3   Cost Effective Lazy Forward [1]

Before, we talk about the algorithmic optimization, I need to introduce the concept of submodular set functions. There are many equivalent definitions of submodular functions, but the one that we will be the most interested in in this post is given below.

Definition 3 (Submodular set function). Let \(\Omega\) be a finite set, then a set function \(f:\mathcal{P}(\Omega)\to\mathbb{R}\) is said to be submodular if for all \(X\subseteq Y\subseteq \Omega\), and for all \(x\notin Y\) \[f(X\cup\{x\})-f(X)\ge f(Y\cup\{x\})-f(Y).\]

Remark. The coverage function \(c\) is a submodular function. And the marginal coverage gain when seen as a function of the set is a monotonic non-increasing function due to the submodularity of \(c\). By monotonic non-increasing, I mean that the values of the marginal gain \(\sigma\) decrease (or stay the same) as we go to larger and larger super-sets as inputs of this function.

Remark. To be more explicit, it is important to note that being a submodular set function, the coverage function must take a set. This is the set \(X\) which is a subset of \(\mathcal{P}(U)\) for some other set \(U\) (see definition 1). Therefore \(X\) is necessarily a set of sets. And therefore, the element \(x\) that we add to \(X\) in definition 3 when we map \(f=c\) is itself a set. The nuance here is that submodular set functions can be functions of any kinds of sets, sets of numbers, sets of cars, sets of python objects, and in this case, sets of sets. And the object \(x\) being added to the set must be of the correct type to make submodularity useful and meaningful.

3.1   The CELF optimization

Observe that all we want is to be able to find the set with the max marginal gain from the leftover set \(L\) in each iteration. As we go on through each iteration, we want to update each set's marginal gain. If the marginal gain's didn't change, we could very easily use a max-heap to get the set with max marginal gain at each iteration. But since \(\sigma\) depends on \(S^t\), it depends on \(t\). Therefore, we need to recompute \(\sigma\) for all sets in each iteration.

However, since \(c\) is submodular, and \(\sigma\) is monotonic non-increasing, we can reduce a lot of recomputations of \(\sigma(S^t,x)\). And we can still use a heap.

Proposition 1. For any iteration \(t\) of the greedy algorithm, \(\sigma^t(x)\le \sigma^{t-1}(x)\).

Proof. This follows directly from \(S^t\supseteq S^{t-1}\), \(\sigma^t(x)=\sigma(S^t,x)\), and \(\sigma^{t-1}(x)=\)\(\sigma(S^{t-1},x)\) and the monotonicity of \(\sigma\).

\(\square\)

Remark. This means \(\sigma\) is monotonous in the usual sense when you see \(\sigma\) as a function of the iteration \(t\).

The following reasoning works for any iteration \(1\lt t\le k\) in algorithm 1, that is we have to have added at least one set to our solution set. We can create a max-heap of all the sets, with the marginal gain as the sorting key. Let us denote this heap by \(h_t\), indexed by the iteration \(t\). Let the max popped elements from these heap be \(h_t^{(0)},h_t^{(1)},\dotsc\) These are just the sets in sorted order of \(\sigma^{t-1}\), that is \(\sigma^{t-1}(h_t^{(i)})\ge \sigma^{t-1}(h_t^{(j)})\) for all \(i\lt j\).

Remark. We won't actually pop them all, but this way it is easier for me to explain. Also, technically, \(h_t^{(i)}\) are sorted by \(\sigma^{\ell}\) where \(\ell\lt t\) because we avoid several recomputations as you will see in the following, but again the above notation makes it easier for me to explain.

Next, in the usual algorithm, we would have computed \(\sigma^{t}\) for all of the sets \(h_t^{(i)}\). However, here we will not need to. We will start computing \(\sigma^{t}\) from the left, that is in the order that we pop, and look for the point \(i\) where \(\sigma^t(h_t^{(i)})\ge\sigma^{t-1}(h_t^{(i+1)})\). Once we find this point, we add the set \(h_t^{(i)}\) to our solution set, and we are done with the iteration. We never compute \(\sigma^t(h_t^{(j)})\) for any \(j\gt i\). In the next iteration, when we must create the heap again, we will use these stale \(\sigma\) values of the leftover sets.

But is it exactly the same solution and if it is then why should this work? This is answered by the following proposition.

Proposition 2. \(h_t^{(i)}\) selected by the above scheme is the same as \(\operatorname{arg\,max}_{x\in L}\sigma^t(x)\).

Proof. Notice that if \(t\gt\ell\), then \(\sigma^t(x)\le\sigma^\ell(x)\). We also know that if \(i\lt j\), then \(\sigma^{t-1}(h_t^{(i)})\ge\sigma^{t-1}(h_t^{(j)})\). Combining the two, we know that the element \(h_t^{(i)}\) such that \(\sigma^t(h_t^{(i)})\ge\sigma^{t-1}(h_t^{(i+1)})\) has the maximum value of \(\sigma^t\).

Why?

1. Because for all \(j\lt i\), we have recomputed \(\sigma^t(h_t^{(j)})\) and know that \(\sigma^t(h_t^{(i)})\gt\sigma^{t-1}(h_t^{(j)})\) (if this wasn't the case, then we would have stopped at an index before \(i\)).
2. We also know that \(\sigma^{t-1}(h_t^{i+1})\ge\sigma^t(h_t^{i+1})\) (due to proposition 1), which can be combined with the stopping condition to get \(\sigma^t(h_t^{(i)})\ge\)\(\sigma^{t-1}(h_t^{(i+1)})\ge\)\(\sigma^{t}(h_t^{(i+1)})\).
3. Finally, for \(j\gt i+1\) we know that, \(\sigma^{t-1}(h_t^{i+1})\ge\sigma^{t-1}(h_t^{j})\) by the way of creation of the heap. Finally, again from proposition 1, we know that \(\sigma^{t-1}(h_t^{(j)})>\sigma^t(h_t^{(j)})\).

Combining the three facts, we get that \(\sigma^t(h_t^{(i)})\ge\sigma^t(h_t^{(j)})\) for all \(j\). Therefore, we still select \(\operatorname{arg\,max}_{x\in L}\sigma^t(x)\), and thus get the exact same solution.

\(\square\)

Remark. The above algorithm works just fine if we replace \(t-1\) with any number \(\lt t\), which is how we represent it in the following algorithm. Further, in the python implementation below, we have used the term stale gain to represent that these values are now any "stale" value of \(\sigma\), and not just \(\sigma^{t-1}\).

Even though we can extract the algorithm from the proof of proposition 2, for the sake of completeness, here is the algorithm.

3.2   Python implementation of the code

This is more of a celf reference for future me, if I ever need a python implementation in future, and forget how to implement this correctly. There are some tricks that I need to pull to make things work since there are some pecularities of python like, for example, python only has min-heap, and there are some added checks to ensure the heap length doesn't become zero after popping an element or that we do not add zero marginal gain elements. Apart from these technical changes, the implementation is almost a verbatim copy of the above algorithm.

import heapq
def max_coverage_with_celf(sets, k=float('inf')):
    covered_elements = set()
    selected_sets = []
    h = [(-len(Set), set_id) for set_id, Set in sets.items()]
    heapq.heapify(h)
    neg_gain, set_id = heapq.heappop(h)
    selected_sets.append(set_id)
    covered_elements.update(sets[set_id])
    ctr = 1
    while h and ctr < k:
        neg_gain, set_id = heapq.heappop(h)
        if neg_gain >= 0:
            break
        new_gain = len(set(sets[set_id]) - covered_elements)
        max_gain, max_set = new_gain, set_id
        if len(h) != 0:
            temp_list = [(new_gain, set_id)]
            stale_gain, peek_h = h[0]
            max_set_idx = 0
            idx = 0
            while new_gain < -stale_gain:
                new_gain = len(set(sets[peek_h]) - covered_elements)
                temp_list.append((new_gain, peek_h))
                idx += 1
                if new_gain > max_gain:
                    max_set_idx = idx
                    max_gain = new_gain
                    max_set = peek_h
                heapq.heappop(h)
                stale_gain, peek_h = h[0]
            temp_list.pop(max_set_idx)
            for new_gain, set_id in temp_list:
                heapq.heappush(h, (-new_gain, set_id))
        selected_sets.append(max_set)
        covered_elements.update(sets[max_set])
        ctr += 1
    return selected_sets




def main():
    S1 = {1, 2, 3}
    S2 = {2, 3, 4}
    S3 = {5, 8}
    S4 = {7, 6}
    sets = {1: S1, 2: S2, 3: S3, 4: S4}
    k = 2
    selected_sets = max_coverage_with_celf(sets, k)
    print(selected_sets)




if __name__=="__main__":
    main()

Concluding remarks

This algorithm is a great example of humans being clever to extract that extra bit of efficiency from the tight-hold of NP-hardness. I like to think that this "trick" is similar to some of those clever tricks that allowed humans to send space-crafts across the solar system as they had to place very precise and energy efficient algorithms in very low resource devices. And efficiency comes from extreme preciseness, you cannot be bothered with any extra computation. However, the speed-up provided by CELF is only meaningful when the distribution of coverages of sets is power-law, or at least very skewed. Finally, I would also like to mention that CELF is not limited to just max coverage, it can also be employed for other problems that have the submodularity structure, like influence maximization, or vertex cover. There might be cases where the function \(\sigma\) might be negative, in which case one must be careful in not selecting those "sets" with "negative" marginal gain, since the objective is maximization.

References

[1] Leskovec, Jure, Andreas Krause, Carlos Guestrin, Christos Faloutsos, Jeanne VanBriesen, and Natalie Glance. "Cost-effective outbreak detection in networks." In Proceedings of the 13th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 420-429. 2007.
[2] Nemhauser, George L., Laurence A. Wolsey, and Marshall L. Fisher. "An analysis of approximations for maximizing submodular set functions—I." Mathematical programming 14 (1978): 265-294.