How can I solve a Set Cover problem in Mathematica

I have $$n$$ non-empty possibly non-disjoint sets $$S_i$$, each having a cost $$c_i$$, and the union $$\Omega=\bigcup_{i=1}^{n}S_i$$. How can I find a selection of $$S_i$$ such that the union is also $$\Omega$$ and the total cost of this selection is minimized?

This is called the unweighted budgeted maximum cover problem.

For example, I have these sets:

(* {set, cost} *)
sets = {
{{1, 2, 6}, 24},
{{3, 5}, 18},
{{1, 5, 6}, 7},
{{4, 5, 6}, 14},
{{2, 3}, 12},
{{1}, 5}
};


The union is Union@@sets[[All,1]] which is {1,2,3,4,5,6}. One could choose the last three sets for a total cost of $$14 + 12 + 5 = 31$$ . Mathematica has FindVertexCover and FindEdgeCover which are graph problems related to the set cover problem. Is there a way to solve set cover problems?

Here's an equivalent version of Daniel's answer that directly uses LinearProgramming:

elements = Union @@ sets[[All, 1]];
res = LinearProgramming[
sets[[All,2]],
ConstantArray[{1,1},Length@elements],
ConstantArray[{0,1},Length@sets],
Integers
]


LinearProgramming::lpip: Warning: integer linear programming will use a machine-precision approximation of the inputs.

{0, 0, 0, 1, 1, 1}

The second argument of the LinearProgramming call is a matrix where the $$n^{\text{th}}$$ row specifies the sets which have $$n$$ as a member:

Transpose @ SparseArray[
] //MatrixForm //TeXForm


$$\begin{array}{cccccc} 1 & 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 \\ \end{array}$$

The third argument specifies that each element must be present at least once, and the fourth argument specifies that either 0 or 1 instances of each of the sets can be used.

The lowest cost is then:

res . sets[[All, 2]]


31

• Thanks. If I solve over the reals I think I can use this for the incremental randomized-rounding solution in case the problem size becomes too large for integer-LP. – flinty Jun 25 at 21:08
• I'm having some trouble with the dimensions when I add another set or change the elements. I tried setting elements = Union @@ sets[[All, 1]]; and then changing the first 6 to Length[sets], the second 6 in the ConstantArray to Length[elements], and the final 6 to Length[sets]. – flinty Jun 25 at 21:29
• @flinty Assuming that the elements range from 1 to max with no missing elements, then the first two 6's should be max, and the last 6 should be Length[sets]. – Carl Woll Jun 26 at 2:02

Could use integer linear programming (ILP). One way to set this up is to use 0-1 "indicator" variables that determine whether a set is used or not in the cover. The objective is to minimize the costs of the used sets. Constraints are that all variables are 0 or 1 and that all elements in the superset appear in sets with an indicator value of 1.

The example can be coded as below.

sets = {{{1, 2, 6}, 24}, {{3, 5}, 18}, {{1, 5, 6}, 7}, {{4, 5, 6},
14}, {{2, 3}, 12}, {{1}, 5}};
vals = Union @@ sets[[All, 1]];
n = Length[vals];
costs = sets[[All, 2]];
vars = Array[c, n];
obj = costs.vars;
c1 = Map[0 <= # <= 1 &, vars];
c2 = Map[Sum[
If[MemberQ[sets[[j, 1]], #], c[j], 0], {j, Length[sets]}] >=
1 &, vals];
constraints = Flatten[{c1, c2}];


Now minimize.

Minimize[{obj, constraints}, vars, Integers]

(* Out[258]= {31, {c[1] -> 0, c[2] -> 0, c[3] -> 0, c[4] -> 1, c[5] -> 1,
c[6] -> 1}} *)


It is straightforward to package this as a function taking a set of sets and costs as argument.

• When the number of elements is different from the number of sets, e.g by adding one more set, the dot product has the wrong dimensions. Please could you suggest a fix for this case. – flinty Jun 25 at 21:35
• My bug. I should be n = Length[sets]; and also use n instead of Length[sets] later on. – Daniel Lichtblau Jun 26 at 15:44

This is how to do it using the greedy algorithm. It keeps selecting sets which give the best ratio of new elements to smallest cost. However, other algorithms may do better. I hope others can suggest improvements or alternative approaches to the problem, such as the LP-relaxation (randomized rounding) method, or perhaps using FindInstance, or even brute-force.

greedyfactor[curunion_, set_] :=
If[ContainsAll[curunion, set[[1]]], Infinity,
set[[2]]/(Length[Union[curunion, set[[1]]]] - Length[curunion])]

greedy[sets_] :=
Module[{target = Union @@ sets[[All, 1]], remaining = sets,
curunion = {}, sel},
Reap[While[! ContainsAll[curunion, target],
sel = First[MinimalBy[remaining, greedyfactor[curunion, #] &]];
Sow[sel];
remaining = DeleteCases[remaining, sel];
curunion = Union[curunion, sel[[1]]]
]][[-1, 1]]
]

sets = {
{{1, 2, 6}, 24},
{{3, 5}, 18},
{{1, 5, 6}, 7},
{{4, 5, 6}, 14},
{{2, 3}, 12},
{{1}, 5}
};

(* confirm we cover all elements *)
chosenSets = greedy[sets]

(* result: {{{1, 5, 6}, 7}, {{2, 3}, 12}, {{4, 5, 6}, 14}} *)

Union @@ chosenSets[[All, 1]] == Union @@ sets[[All, 1]]

(* result: True *)

(* get the final cost *)
cost = Total[chosenSets[[All, 2]]]

(* result: 33 *)