MLIR — A Global Optimization and Dataflow Analysis

Table of Contents

In this article we’ll implement a global optimization pass, and show how to use the dataflow analysis framework to verify the results of our optimization.

The code for this article is in this pull request, and as usual the commits are organized to be read in order.

The noisy arithmetic problem

This demonstration is based on a simplified model of computation relevant to the HEIR project. You don’t need to be familiar with that project to follow this article, but if you’re wondering why someone would ever want the kind of optimization I’m going to write, that project is why.

The basic model is “noisy integer arithmetic.” That is, a program can have integer types of bounded width, and each integer is associated with some unknown (but bounded with high probability) symmetric random noise. You can imagine the “real” integer being the top 5 bits of a 32-bit integer, and the bottom 27 bits storing the noise. When a new integer is created, it magically has a random signed 12-bit integer added to it. When you apply operations to combine integers, the noise grows. Adding two integers adds their noise, and at worst you get one more bit of noise. Scaling an integer by a statically-known constant scales the noise by a constant. Multiplying two integers multiplies their noise values, and you get twice the bits of noise. As long as your program stays below 27 bits of noise, you can still “recover” the original 5-bit integer at the end of the program. Such a program is called legal, and otherwise, the output is random junk and the program is called illegal.

Finally, there is an expensive operation called reduce_noise that can explicitly reduce the noise of a noisy integer back to the base level of 12 bits. This operation has a statically known cost relative to the standard operations.

Note that starting from two noisy integers, each with 12 bits of noise, a single multiplication op brings you jarringly close to ruin. You would have at most 24 bits of noise, which is close to the maximum of 26. But the input IR we start with may have arbitrary computations that do not respect the noise limits. The goal of our optimization is to rewrite an MLIR region so that the noisy integer math never exceeds the noise limit at any given step.

A trivial way to do that would be to insert reduce_noise ops greedily, whenever an op would bring the noise of a value too high. Inserting such ops may be necessary, but a more suitable goal would be to minimize the overall cost of the program, subject to the constraint that the noisy arithmetic is legal. One could do this by inserting reduce_noise ops more strategically, or by rewriting the program to reduce the need for reduce_noise ops, or both. We’ll focus on the former: finding the best place to insert reduce_noise ops without rewriting the rest of the program.

The noisy dialect

We previously wrote about defining a new dialect, and the noisy dialect we created for this article has little new to show. This commit defines the types and ops, hard coding the 32-bit width and 5-bit semantic input type for simplicity, as well as the values of 12 bits of initial noise and 26 bits of max noise.

Note that the noise bound is not expressed on the type as an attribute. If it were, we’d run into a few problems: first, whenever you insert a reduce_noise op, you’d have to update the types on all the downstream ops. Second, it would prevent you from expressing control flow, since the noise bound cannot be statically inferred from the source code when there are two possible paths that could result in different noise values.

So instead, we need a way to compute the noise values, and associate them with each SSA value, and account for control flow. This is what an analysis pass is designed to do.

An analysis pass is just a class

The typical use of an analysis pass is to construct a data structure that encodes global information about a program, which can then be re-used during different parts of a pass. I imagined there would be more infrastructure around analysis passes in MLIR, but it’s quite simple. You define a C++ class with a constructor that takes an Operation *, and construct it basically whenever you need it. The only infrastructure for it involves storing and caching the constructed analysis within a pass, and understanding when an analysis needs to be recomputed (always between passes, by default).

By way of example, over at the HEIR project I made a simple analysis that chooses a unique variable name for every SSA value in a program, which I then used to generate code in an output language that needed variable names.

For this article we’ll see two analysis passes. One will formulate and solve the optimization problem that decides where to insert reduce_noise operations. This will be one of the “class that does anything” kind of analysis pass. The other analysis pass will rely on MLIR’s data flow analysis framework to propagate the noise model through the IR. This one will actually not require us to write an analysis from scratch, but instead will be implemented by means of the existing IntegerRangeAnalysis, which only requires us to implement an interface on each op that describes how the op affects the noise. This will be used in our pass to verify that the inserted reduce_noise operations ensure, if nothing else, that the noise never exceeds the maximum allowable noise.

We’ll start with the data flow analysis.

Reusing IntegerRangeAnalysis

Data flow analysis is a classical static analysis technique for propagating information through a program’s IR. It is one part of Frances Allen’s Turing Award. This article gives a good introduction and additional details, but I will paraphrase it briefly here in the context of IntegerRangeAnalysis.

The basic idea is that you want to get information about what possible values an integer-typed value can have at any point in a given program. If you see x = 7, then you know exactly what x is. If you see something like

func (%x : i8) { 
  %1 = arith.extsi %x : i8 to i32 
  %2 = arith.addi %x, %x : i32

then you know that %2 can be at most a signed 9-bit integer, because it started as an 8-bit integer, and adding two such integers together can’t fill up more than one extra bit.

In such cases, one can find optimizations, like the int-range-optimizations pass in MLIR, which looks at comparison ops arith.cmpi and determines if it can replace them with constants. It does this by looking at the integer range analysis for the two operands. E.g., given the op x > y, if you know y‘s maximum value is less than x‘s minimum value, then you can replace it with a constant true.

Computing the data flow analysis requires two ingredients called a transfer function and a join operation. The transfer function describes what the output integer range should be for a given op and a given set of input integer ranges. This can be an arbitrary function. The join operation describes how to combine two or more integer ranges when you get to a point at the program in which different branches of control flow merge. For example,

def fn(branch):
  x = 7
  if branch:
     y = x * x
     y = 2*x
  return y

The value of y just before returning cannot be known exactly, but in one branch you know it’s 14, and in another it’s 49. So the final value of y could be estimated as being in the range [14, 49]. Here the join function computes the smallest integer range containing both estimates. [Aside: it could instead use the set {14, 49} to be more precise, but that is not what IntegerRangeAnalysis happens to do]

In order for a data flow analysis to work properly, the values being propagated and the join function must together form a semilattice, which is a partially-ordered set in which every two elements have a least upper bound, that upper bound is computed by join, and join itself must be associative, commutative, and idempotent. For dataflow analysis, the semilattice must also be finite. This is often expressed by having distinct “top” and “bottom” elements as defaults. “Top” represents “could be anything,” and sometimes expresses that a more precise bound would be too computationally expensive to continue to track. “Bottom” usually represents an uninitialized value.

Once you have this, then MLIR provides a general algorithm to propagate values through the IR via a technique called Kildall’s method, which iteratively updates the SSA values, applying the transfer function and joining at the merging of control flow paths, until the process reaches a fixed point.

Here are MLIR’s official docs on dataflow analysis, and here is the RFC where the current data flow solver framework was introduced. In our situation, we want to use the solver framework with the existing IntegerRangeAnalysis, which only asks that we implement the transfer function by implementing InferIntRangeInterface on our ops.

This commit does just that. This requires adding the DeclareOpInterfaceMethods<InferIntRangeInterface> to all relevant ops. That in turn generates function declarations for

void MyOp::inferResultRanges(
    ArrayRef<ConstantIntRanges> inputRanges, SetIntRangeFn setResultRange);

The ConstantIntRange is a dataclass holding a min and max integer value. inputRanges represents the known bounds on the inputs to the operation in question, and SetIntRangeFn is the callback used to produce the result.

For example, for AddOp we can implement it as

ConstantIntRanges unionPlusOne(ArrayRef<ConstantIntRanges> inputRanges) {
  auto lhsRange = inputRanges[0];
  auto rhsRange = inputRanges[1];
  auto joined = lhsRange.rangeUnion(rhsRange);
  return ConstantIntRanges::fromUnsigned(joined.umin(), joined.umax() + 1);

void AddOp::inferResultRanges(ArrayRef<ConstantIntRanges> inputRanges,
                              SetIntRangeFn setResultRange) {
  setResultRange(getResult(), unionPlusOne(inputRanges));

A MulOp is similarly implemented by summing the maxes. Meanwhile, EncodeOp and ReduceNoiseOp each set the initial range to [0, 12]. So the min will always be zero, and we really only care about the max.

The next commit defines an empty pass that will contain our analyses and optimizations, and this commit shows how the integer range analysis is used to validate an IR’s noise growth. In short, you load the IntegerRangeAnalysis and its dependent DeadCodeAnalysis, run the solver, and then walk the IR, asking the solver via lookupState to give the resulting value range for each op’s result, and comparing it against the maximum.

void runOnOperation() {
    Operation *module = getOperation();

    DataFlowSolver solver;
    if (failed(solver.initializeAndRun(module)))

    auto result = module->walk([&](Operation *op) {
      if (!llvm::isa<noisy::AddOp, noisy::SubOp, noisy::MulOp,
                     noisy::ReduceNoiseOp>(*op)) {
        return WalkResult::advance();
      const dataflow::IntegerValueRangeLattice *opRange =
      if (!opRange || opRange->getValue().isUninitialized()) {
            << "Found op without a set integer range; did the analysis fail?";
        return WalkResult::interrupt();

      ConstantIntRanges range = opRange->getValue().getValue();
      if (range.umax().getZExtValue() > MAX_NOISE) {
        op->emitOpError() << "Found op after which the noise exceeds the "
                             "allowable maximum of "
                          << MAX_NOISE
                          << "; it was: " << range.umax().getZExtValue()
                          << "\n";
        return WalkResult::interrupt();

      return WalkResult::advance();

    if (result.wasInterrupted())

Finally, in this commit we add a test that exercises it:

func.func @test_op_syntax() -> i5 {
  %0 = arith.constant 3 : i5
  %1 = arith.constant 4 : i5
  %2 = noisy.encode %0 : i5 -> !noisy.i32
  %3 = noisy.encode %1 : i5 -> !noisy.i32
  %4 = noisy.mul %2, %3 : !noisy.i32
  %5 = noisy.mul %4, %4 : !noisy.i32
  %6 = noisy.mul %5, %5 : !noisy.i32
  %7 = noisy.mul %6, %6 : !noisy.i32
  %8 = noisy.decode %7 : !noisy.i32 -> i5
  return %8 : i5

Running tutorial-opt --noisy-reduce-noise on this file produces the following error:

 error: 'noisy.mul' op Found op after which the noise exceeds the allowable maximum of 26; it was: 48

  %5 = noisy.mul %4, %4 : !noisy.i32
mlir-tutorial/tests/noisy_reduce_noise.mlir:11:8: note: see current operation: %5 = "noisy.mul"(%4, %4) : (!noisy.i32, !noisy.i32) -> !noisy.i32

And if you run in debug mode with --debug --debug-only=int-range-analysis, you will see the per-op propagations printed to the terminal

$  bazel run tools:tutorial-opt -- --noisy-reduce-noise-optimizer $PWD/tests/noisy_reduce_noise.mlir --debug --debug-only=int-range-analysis
Inferring ranges for %c3_i5 = arith.constant 3 : i5
Inferred range unsigned : [3, 3] signed : [3, 3]
Inferring ranges for %c4_i5 = arith.constant 4 : i5
Inferred range unsigned : [4, 4] signed : [4, 4]
Inferring ranges for %0 = noisy.encode %c3_i5 : i5 -> !noisy.i32
Inferred range unsigned : [0, 12] signed : [0, 12]
Inferring ranges for %1 = noisy.encode %c4_i5 : i5 -> !noisy.i32
Inferred range unsigned : [0, 12] signed : [0, 12]
Inferring ranges for %2 = noisy.mul %0, %1 : !noisy.i32
Inferred range unsigned : [0, 24] signed : [0, 24]
Inferring ranges for %3 = noisy.mul %2, %2 : !noisy.i32
Inferred range unsigned : [0, 48] signed : [0, 48]
Inferring ranges for %4 = noisy.mul %3, %3 : !noisy.i32
Inferred range unsigned : [0, 96] signed : [0, 96]
Inferring ranges for %5 = noisy.mul %4, %4 : !noisy.i32
Inferred range unsigned : [0, 192] signed : [0, 192]

As a quick aside, there was one minor upstream problem preventing me from reusing IntegerRangeAnalysis, which I patched in This means I also had to update the LLVM commit hash used by this project in this commit.

An ILP optimization pass

Next, we build an analysis that solves a global optimization to insert reduce_noise ops efficiently. As mentioned earlier, this is a “do anything” kind of analysis, so we put all of the logic into the analysis’s construtor.

[Aside: I wouldn’t normally do this, because constructors don’t have return values so it’s hard to signal failure; but the API for the analysis specifies the constructor takes as input the Operation * to analyze, and I would expect any properly constructed object to be “ready to use.” Maybe someone who knows C++ better will comment and shed some wisdom for me.]

This commit sets up the analysis shell and interface.

class ReduceNoiseAnalysis {
  ReduceNoiseAnalysis(Operation *op);
  ~ReduceNoiseAnalysis() = default;

  /// Return true if a reduce_noise op should be inserted after the given
  /// operation, according to the solution to the optimization problem.
  bool shouldInsertReduceNoise(Operation *op) const {
    return solution.lookup(op);

  llvm::DenseMap<Operation *, bool> solution;

This commit adds a workspace dependency on Google’s or-tools package (“OR” stands for Operations Research here, a.k.a. discrete optimization), which comes bundled with a number of nice solvers, and an API for formulating optimization problems. And this commit implements the actual solver model.

Now this model is quite a bit of code, and this article is not the best place to give a full-fledged introduction to linear programming, modeling techniques, or the OR-tools API. What I’ll do instead is explain the model in detail here, give a few small notes on how that translates to the OR-tools C++ API. If you want a gentler background on linear programming, see my article series about diet optimization (part 1, part 2).

All linear programs specify a linear function as an objective to minimize, along with a set of linear equalities and inequalities that constrain the solution. In a standard linear program, the variables must be continuously valued. In a mixed-integer linear program, some of those variables are allowed to be discrete integers, which, it turns out, makes it possible to solve many more problems, but requires completely different optimization techniques and may result in exponentially slow runtime. So many techniques in operations research relate to modeling a problem in such a way that the number of integer variables is relatively small.

Our linear model starts by defining some basic variables. Some variables in the model represent “decisions” that we can make, and others represent “state” that reacts to the decisions via constraints.

  • For each operation $x$, a $\{0, 1\}$-valued variable $\textup{InsertReduceNoise}_x$. Such a variable is 1 if and only if we insert a reduce_noise op after the operation $x$.
  • For each SSA-value $v$ that is input or output to a noisy op, a continuous-valued variable $\textup{NoiseAt}_v$. This represents the upper bound of the noise at value $v$.

In particular, the solver’s performance will get worse as the number of binary variables increases, which in this case corresponds to the number of noisy ops.

The objective function, with a small caveat explained later, is simply the sum of the decision variables, and we’d like to minimize it. Each reduce_noise op is considered equally expensive, and there is no special nuance here about scheduling them in parallel or in serial.

Next, we add constraints. First, $0 \leq \textup{NoiseAt}_v \leq 26$, which asserts that no SSA value can exceed the max noise. Second, we need to enforce that an encode op fixes the noise of its output to 12, i.e., for each encode op $x$, we add the constraint $\textup{NoiseAt}_{\textup{result}(x)} = 12$.

Finally, we need constraints that say that if you choose to insert a reduce_noise op, then the noise is reset to 12, otherwise it is set to the appropriate function of the inputs. This is where the modeling gets a little tricky, but multiplication is easier so let’s start there.

Fix a multiplication op $x$, its two input SSA values $\textup{LHS}, \textup{RHS}$, and its output $\textup{RES}$. As a piecewise function, we want a constraint like:

\[ \textup{NoiseAt}_\textup{RES} = \begin{cases} \textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS} & \text{ if } \textup{InsertReduceNoise}_x = 0 \\ 12 & \text{ if } \textup{InsertReduceNoise}_x = 1 \end{cases} \]

This isn’t linear, but we can combine the two branches to

\[ \begin{aligned} \textup{NoiseAt}_\textup{RES} &= (1 – \textup{ InsertReduceNoise}_x) (\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS}) \\ & + 12 \textup{ InsertReduceNoise}_x \end{aligned} \]

This does the classic trick of using a bit as a controlled multiplexer, but it’s still not linear. We can make it linear, however, by replacing this one constraint with four constraints, and an auxiliary constant $C=100$ that we know is larger than the possible range of values that the $\textup{NoiseAt}_v$ variables can attain. Those four linear constraints are:

\[ \begin{aligned}
\textup{NoiseAt}_\textup{RES} &\geq 12 \textup{ InsertReduceNoise}_x \\
\textup{NoiseAt}_\textup{RES} &\leq 12 + C(1 – \textup{InsertReduceNoise}_x) \\
\textup{NoiseAt}_\textup{RES} &\geq (\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS}) – C \textup{ InsertReduceNoise}_x \\
\textup{NoiseAt}_\textup{RES} &\leq (\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS}) + C \textup{ InsertReduceNoise}_x \\
\end{aligned} \]

Setting the decision variable to zero results in the first two equations being trivially satisfied. Setting it to 1 causes the first two equations to be equivalent to $\textup{NoiseAt}_\textup{RES} = 12$. Likewise, the second two constraints are trivial when the decision variable is 1, and force the output noise to be equal to the sum of the two input noises when set to zero.

The addition op is handled similarly, except that the term $(\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS})$ is replaced by something non-linear, namely $1 + \max(\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS})$. We can still handle that, but it requires an extra modeling trick. We introduce a new variable $Z_x$ for each add op $x$, and two constraints:

\[ \begin{aligned} Z_v &\geq 1 + \textup{NoiseAt}_{LHS} \\ Z_v &\geq 1 + \textup{NoiseAt}_{RHS} \end{aligned} \]

Together these ensure that $Z_v$ is at least 1 plus the max of the two input noises, but it doesn’t force equality. To achieve that, we add $Z_v$ to the minimization objective (alongside the sum of the decision variables) with a small penalty to ensure the solver tries to minimize them. Since they have trivially minimal values equal to “1 plus the max,” the solver will have no trouble optimizing them, and this will be effectively an equality constraint.

[Aside: Whenever you do this trick, you have to convince yourself that the solver won’t somehow be able to increase $Z_v$ as a trade-off against lower values of other objective terms, and produce a lower overall objective value. Solvers are mischievous and cannot be trusted. In our case, there is no risk: if you were to increase $Z_v$ below its minimum value, that would only increase the noise propagation through add ops, meaning the solver would have to compensate by potentially adding even more reduce_noise ops!]

Then, the constraint for an add op uses $Z_v$ in place of $(\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS})$ the mul op.

The only other minor aspect of this solver model is that these constraints enforce the consistency of the noise propagation after a reduce_noise op may be added, but if a reduce_noise op is added, it doesn’t necessarily enforce the noise growth of the output of the op before it’s input to reduce_noise. We can achieve this by adding new constraints expressing $(\textup{NoiseAt}_{LHS} + \textup{NoiseAt}_{RHS}) leq 26$ and $Z_v \leq 26$ for multiplication and addition ops, respectively.

When converting this to the OR-tools C++ API, as we did in this commit, a few minor things to note:

  • You can specify upper and lower bounds on a variable at variable creation time, rather than as separate constraints. You’ll see this in solver->MakeNumVar(min, max, name).
  • Constraints must be specified in the form min <= expr <= max, where min and max are constants and expr is a linear combination of variables, meaning that one has to manually re-arrange and simplify all the equations above so the variables are all on one side and the constants on the other. (The OR-tools Python API is more expressive, but we don’t have it here.)
  • The constraints and the objective are specified by SetCoefficient, which sets the coefficient of a variable in a linear combination one at a time.

Finally, this commit implements the part of the pass that uses the solver’s output to insert new reduce_noise ops. And this commit adds some more tests.

An example of its use:

// This test checks that the solver can find a single insertion point
// for a reduce_noise op that handles two branches, each of which would
// also need a reduce_noise op if handled separately.
func.func @test_single_insertion_branching() -> i5 {
  %0 = arith.constant 3 : i5
  %1 = arith.constant 4 : i5
  %2 = noisy.encode %0 : i5 -> !noisy.i32
  %3 = noisy.encode %1 : i5 -> !noisy.i32
  // Noise: 12
  %4 = noisy.mul %2, %3 : !noisy.i32
  // Noise: 24

  // branch 1
  %b1 = noisy.add %4, %3 : !noisy.i32
  // Noise: 25
  %b2 = noisy.add %b1, %3 : !noisy.i32
  // Noise: 25
  %b3 = noisy.add %b2, %3 : !noisy.i32
  // Noise: 26
  %b4 = noisy.add %b3, %3 : !noisy.i32
  // Noise: 27

  // branch 2
  %c1 = noisy.sub %4, %2 : !noisy.i32
  // Noise: 25
  %c2 = noisy.sub %c1, %3 : !noisy.i32
  // Noise: 25
  %c3 = noisy.sub %c2, %3 : !noisy.i32
  // Noise: 26
  %c4 = noisy.sub %c3, %3 : !noisy.i32
  // Noise: 27

  %x1 = noisy.decode %b4 : !noisy.i32 -> i5
  %x2 = noisy.decode %c4 : !noisy.i32 -> i5
  %x3 = arith.addi %x1, %x2 : i5
  return %x3 : i5

And the output:

  func.func @test_single_insertion_branching() -> i5 {
    %c3_i5 = arith.constant 3 : i5
    %c4_i5 = arith.constant 4 : i5
    %0 = noisy.encode %c3_i5 : i5 -> !noisy.i32
    %1 = noisy.encode %c4_i5 : i5 -> !noisy.i32
    %2 = noisy.mul %0, %1 : !noisy.i32
    %3 = noisy.reduce_noise %2 : !noisy.i32
    %4 = noisy.add %3, %1 : !noisy.i32
    %5 = noisy.add %4, %1 : !noisy.i32
    %6 = noisy.add %5, %1 : !noisy.i32
    %7 = noisy.add %6, %1 : !noisy.i32
    %8 = noisy.sub %3, %0 : !noisy.i32
    %9 = noisy.sub %8, %1 : !noisy.i32
    %10 = noisy.sub %9, %1 : !noisy.i32
    %11 = noisy.sub %10, %1 : !noisy.i32
    %12 = noisy.decode %7 : !noisy.i32 -> i5
    %13 = noisy.decode %11 : !noisy.i32 -> i5
    %14 = arith.addi %12, %13 : i5
    return %14 : i5

Optimization Models for Subset Cover

In a recent newsletter article I complained about how researchers mislead about the applicability of their work. I gave SAT solvers as an example. People provided interesting examples in response, but what was new to me was the concept of SMT (Satisfiability Modulo Theories), an extension to SAT. SMT seems to have more practical uses than vanilla SAT (see the newsletter for details).

I wanted to take some time to explore SMT solvers, and I landed on Z3, an open-source SMT solver from Microsoft. In particular, I wanted to compare it to ILP (Integer Linear Programing) solvers, which I know relatively well. I picked a problem that I thought would work better for SAT-ish solvers than for ILPs: subset covering (explained in the next section). If ILP still wins against Z3, then that would be not so great for the claim that SMT is a production strength solver.

All the code used for this post is on Github.

Subset covering

A subset covering is a kind of combinatorial design, which can be explained in terms of magic rings.

An adventurer stumbles upon a chest full of magic rings. Each ring has a magical property, but some pairs of rings, when worn together on the same hand, produce a combined special magical effect distinct to that pair.

The adventurer would like to try all pairs of rings to catalogue the magical interactions. With only five fingers, how can we minimize the time spent trying on rings?

Mathematically, the rings can be described as a set $ X$ of size $ n$. We want to choose a family $ F$ of subsets of $ X$, with each subset having size 5 (five fingers), such that each subset of $ X$ of size 2 (pairs of rings) is contained in some subset of $ F$. And we want $ F$ to be as small as possible.

Subset covering is not a “production worthy” problem. Rather, I could imagine it’s useful in some production settings, but I haven’t heard of one where it is actually used. I can imagine, for instance, that a cluster of machines has some bug occurring seemingly at random for some point-to-point RPCs, and in tracking down the problem, you want to deploy a test change to subsets of servers to observe the bug occurring. Something like an experiment design problem.

If you generalize the “5” in “5 fingers” to an arbitrary positive integer $ k$, and the “2” in “2 rings” to $ l < k$, then we have the general subset covering problem. Define $ M(n, k, l)$ to be the minimal number of subsets of size $ k$ needed to cover all subsets of size $ l$. This problem was studied by Erdős, with a conjecture subsequently proved by Vojtěch Rödl, that asymptotically $ M(n,k,l)$ grows like $ \binom{n}{l} / \binom{k}{l}$. Additional work by Joel Spencer showed that a greedy algorithm is essentially optimal.

However, all of the constructive algorithms in these proofs involve enumerating all $ \binom{n}{k}$ subsets of $ X$. This wouldn’t scale very well. You can alternatively try a “random” method, incurring a typically $ \log(r)$ factor of additional sets required to cover a $ 1 – 1/r$ fraction of the needed subsets. This is practical, but imperfect.

To the best of my knowledge, there is no exact algorithm, that both achieves the minimum and is efficient in avoiding constructing all $ \binom{n}{k}$ subsets. So let’s try using an SMT solver. I’ll be using the Python library for Z3.

Baseline: brute force Z3

For a baseline, let’s start with a simple Z3 model that enumerates all the possible subsets that could be chosen. This leads to an exceedingly simple model to compare the complex models against.

Define boolean variables $ \textup{Choice}_S$ which is true if and only if the subset $ S$ is chosen (I call this a “choice set”). Define boolean variables $ \textup{Hit}_T$ which is true if the subset $ T$ (I call this a “hit set”) is contained in a chosen choice set. Then the subset cover problem can be defined by two sets of implications.

First, if $ \textup{Choice}_S$ is true, then so must all $ \textup{Hit}_T$ for $ T \subset S$. E.g., for $ S = \{ 1, 2, 3 \}$ and $ l=2$, we get

$ \displaystyle \begin{aligned} \textup{Choice}_{(1,2,3)} &\implies \textup{Hit}_{(1,2)} \\ \textup{Choice}_{(1,2,3)} &\implies \textup{Hit}_{(1,3)} \\ \textup{Choice}_{(1,2,3)} &\implies \textup{Hit}_{(2,3)} \end{aligned}$

In Python this looks like the following (note this program has some previously created lookups and data structures containing the variables)

for choice_set in choice_sets:
  for hit_set_key in combinations(choice_set.elements, hit_set_size):
    hit_set = hit_set_lookup[hit_set_key]
      z3.Implies(choice_set.variable, hit_set.variable))

Second, if $ \textup{Hit}_T$ is true, it must be that some $ \textup{Choice}_S$ is true for some $ S$ containing $ T$ as a subset. For example,

$ \displaystyle \begin{aligned} \textup{Hit}_{(1,2)} &\implies \\ & \textup{Choice}_{(1,2,3,4)} \textup{ OR} \\ & \textup{Choice}_{(1,2,3,5)} \textup{ OR} \\ & \textup{Choice}_{(1,2,4,5)} \textup{ OR } \cdots \\ \end{aligned}$

In code,

for hit_set in hit_sets.values():
  relevant_choice_set_vars = [
    for choice_set in hit_set_to_choice_set_lookup[hit_set]

Next, in this experiment we’re allowing the caller to specify the number of choice sets to try, and the solver should either return SAT or UNSAT. From that, we can use a binary search to find the optimal number of sets to pick. Thus, we have to limit the number of $ \textup{Choice}_S$ that are allowed to be true and false. Z3 supports boolean cardinality constraints, apparently with a special solver to handle problems that have them. Otherwise, the process of encoding cardinality constraints as SAT formulas is not trivial (and the subject of active research). But the code is simple enough:

args = [cs.variable for cs in choice_sets] + [parameters.num_choice_sets]
choice_sets_at_most = z3.AtMost(*args)
choice_sets_at_least = z3.AtLeast(*args)

Finally, we must assert that every $ \textup{Hit}_T$ is true.

solver = z3.Solver()
for hit_set in hit_sets.values():

for impl in implications:


Running it for $ n=7, k=3, l=2$, and seven choice sets (which is optimal), we get

>>> SubsetCoverZ3BruteForce().solve(
[(0, 1, 3), (0, 2, 4), (0, 5, 6), (1, 2, 6), (1, 4, 5), (2, 3, 5), (3, 4, 6)]
SubsetCoverSolution(status=<SolveStatus.SOLVED: 1>, solve_time_seconds=0.018305063247680664)

Interestingly, Z3 refuses to solve marginally larger instances. For instance, I tried the following and Z3 times out around $ n=12, k=6$ (about 8k choice sets):

from math import comb

for n in range(8, 16):
  k = int(n / 2)
  l = 3
  max_num_sets = int(2 * comb(n, l) / comb(k, l))
  params = SubsetCoverParameters(


After taking a long time to generate the larger models, Z3 exceeds my 15 minute time limit, suggesting exponential growth:

status               solve_time_seconds  num_elements  choice_set_size  hit_set_size  num_choice_sets
SolveStatus.SOLVED   0.0271              8             4                3             28
SolveStatus.SOLVED   0.0346              9             4                3             42
SolveStatus.SOLVED   0.0735              10            5                3             24
SolveStatus.SOLVED   0.1725              11            5                3             33
SolveStatus.SOLVED   386.7376            12            6                3             22
SolveStatus.UNKNOWN  900.1419            13            6                3             28
SolveStatus.UNKNOWN  900.0160            14            7                3             20
SolveStatus.UNKNOWN  900.0794            15            7                3             26

An ILP model

Next we’ll see an ILP model for the sample problem. Note there are two reasons I expect the ILP model to fall short. First, the best solver I have access to is SCIP, which, despite being quite good is, in my experience, about an order of magnitude slower than commercial alternatives like Gurobi. Second, I think this sort of problem seems to not be very well suited to ILPs. It would take quite a bit longer to explain why (maybe another post, if you’re interested), but in short well-formed ILPs have easily found feasible solutions (this one does not), and the LP-relaxation of the problem should be as tight as possible. I don’t think my formulation is very tight, but it’s possible there is a better formulation.

Anyway, the primary difference in my ILP model from brute force is that the number of choice sets is fixed in advance, and the members of the choice sets are model variables. This allows us to avoid enumerating all choice sets in the model.

In particular, $ \textup{Member}_{S,i} \in \{ 0, 1 \}$ is a binary variable that is 1 if and only if element $ i$ is assigned to be in set $ S$. And $ \textup{IsHit}_{T, S} \in \{0, 1\}$ is 1 if and only if the hit set $ T$ is a subset of $ S$. Here “$ S$” is an index over the subsets, rather than the set itself, because we don’t know what elements are in $ S$ while building the model.

For the constraints, each choice set $ S$ must have size $ k$:

$ \displaystyle \sum_{i \in X} \textup{Member}_{S, i} = k$

Each hit set $ T$ must be hit by at least one choice set:

$ \displaystyle \sum_{S} \textup{IsHit}_{T, S} \geq 1$

Now the tricky constraint. If a hit set $ T$ is hit by a specific choice set $ S$ (i.e., $ \textup{IsHit}_{T, S} = 1$) then all the elements in $ T$ must also be members of $ S$.

$ \displaystyle \sum_{i \in T} \textup{Member}_{S, i} \geq l \cdot \textup{IsHit}_{T, S}$

This one works by the fact that the left-hand side (LHS) is bounded from below by 0 and bounded from above by $ l = |T|$. Then $ \textup{IsHit}_{T, S}$ acts as a switch. If it is 0, then the constraint is vacuous since the LHS is always non-negative. If $ \textup{IsHit}_{T, S} = 1$, then the right-hand side (RHS) is $ l = |T|$ and this forces all variables on the LHS to be 1 to achieve it.

Because we fixed the number of choice sets as a parameter, the objective is 1, and all we’re doing is looking for a feasible solution. The full code is here.

On the same simple example as the brute force

>>> SubsetCoverILP().solve(
[(0, 1, 3), (0, 2, 6), (0, 4, 5), (1, 2, 4), (1, 5, 6), (2, 3, 5), (3, 4, 6)]
SubsetCoverSolution(status=<SolveStatus.SOLVED: 1>, solve_time_seconds=0.1065816879272461)

It finds the same solution in about 10x the runtime as the brute force Z3 model, though still well under one second.

On the “scaling” example, it fares much worse. With a timeout of 15 minutes, it solves $ n=8,$ decently fast, $ n=9,12$ slowly, and times out on the rest.

status               solve_time_seconds  num_elements  choice_set_size  hit_set_size  num_choice_sets
SolveStatus.SOLVED   1.9969              8             4                3             28
SolveStatus.SOLVED   306.4089            9             4                3             42
SolveStatus.UNKNOWN  899.8842            10            5                3             24
SolveStatus.UNKNOWN  899.4849            11            5                3             33
SolveStatus.SOLVED   406.9502            12            6                3             22
SolveStatus.UNKNOWN  902.7807            13            6                3             28
SolveStatus.UNKNOWN  900.0826            14            7                3             20
SolveStatus.UNKNOWN  900.0731            15            7                3             26

A Z3 Boolean Cardinality Model

The next model uses Z3. It keeps the concept of Member and Hit variables, but they are boolean instead of integer. It also replaces the linear constraints with implications. The constraint that forces a Hit set’s variable to be true when some Choice set contains all its elements is (for each $ S, T$)

$ \displaystyle \left ( \bigwedge_{i \in T} \textup{Member}_{S, i} \right ) \implies \textup{IsHit}_T$

Conversely, A Hit set’s variable being true implies its members are in some choice set.

$ \displaystyle \textup{IsHit}_T \implies \bigvee_{S} \bigwedge_{i \in T} \textup{Member}_{S, i}$

Finally, we again use boolean cardinality constraints AtMost and AtLeast so that each choice set has the right size.

The results are much better than the ILP: it solves all of the instances in under 3 seconds

status              solve_time_seconds  num_elements  choice_set_size  hit_set_size  num_choice_sets
SolveStatus.SOLVED  0.0874              8             4                3             28
SolveStatus.SOLVED  0.1861              9             4                3             42
SolveStatus.SOLVED  0.1393              10            5                3             24
SolveStatus.SOLVED  0.2845              11            5                3             33
SolveStatus.SOLVED  0.2032              12            6                3             22
SolveStatus.SOLVED  1.3661              13            6                3             28
SolveStatus.SOLVED  0.8639              14            7                3             20
SolveStatus.SOLVED  2.4877              15            7                3             26

A Z3 integer model

Z3 supports implications on integer equation equalities, so we can try a model that leverages this by essentially converting the boolean model to one where the variables are 0-1 integers, and the constraints are implications on equality of integer formulas (all of the form “variable = 1”).

I expect this to perform worse than the boolean model, even though the formulation is almost identical. The details of the model are here, and it’s so similar to the boolean model above that it needs no extra explanation.

The runtime is much worse, but surprisingly it still does better than the ILP model.

status              solve_time_seconds  num_elements  choice_set_size  hit_set_size  num_choice_sets
SolveStatus.SOLVED  2.1129              8             4                3             28
SolveStatus.SOLVED  14.8728             9             4                3             42
SolveStatus.SOLVED  7.6247              10            5                3             24
SolveStatus.SOLVED  25.0607             11            5                3             33
SolveStatus.SOLVED  30.5626             12            6                3             22
SolveStatus.SOLVED  63.2780             13            6                3             28
SolveStatus.SOLVED  57.0777             14            7                3             20
SolveStatus.SOLVED  394.5060            15            7                3             26

Harder instances

So far all the instances we’ve been giving the solvers are “easy” in a sense. In particular, we’ve guaranteed there’s a feasible solution, and it’s easy to find. We’re giving roughly twice as many sets as are needed. There are two ways to make this problem harder. One is to test on unsatisfiable instances, which can be harder because the solver has to prove it can’t work. Another is to test on satisfiable instances that are hard to find, such as those satisfiable instances where the true optimal number of choice sets is given as the input parameter. The hardest unsatisfiable instances are also the ones where the number of choice sets allowed is one less than optimal.

Let’s test those situations. Since $ M(7, 3, 2) = 7$, we can try with 7 choice sets and 6 choice sets.

For 7 choice sets (the optimal value), all the solvers do relatively well

method                    status  solve_time_seconds  num_elements  choice_set_size  hit_set_size  num_choice_sets
SubsetCoverILP            SOLVED  0.0843              7             3                2             7
SubsetCoverZ3Integer      SOLVED  0.0938              7             3                2             7
SubsetCoverZ3BruteForce   SOLVED  0.0197              7             3                2             7
SubsetCoverZ3Cardinality  SOLVED  0.0208              7             3                2             7

For 6, the ILP struggles to prove it’s infeasible, and the others do comparatively much better (at least 17x better).

method                    status      solve_time_seconds  num_elements  choice_set_size  hit_set_size  num_choice_sets
SubsetCoverILP            INFEASIBLE  120.8593            7             3                2             6
SubsetCoverZ3Integer      INFEASIBLE  3.0792              7             3                2             6
SubsetCoverZ3BruteForce   INFEASIBLE  0.3384              7             3                2             6
SubsetCoverZ3Cardinality  INFEASIBLE  7.5781              7             3                2             6

This seems like hard evidence that Z3 is better than ILPs for this problem (and it is), but note that the same test on $ n=8$ fails for all models. They can all quickly prove $ 8 < M(8, 3, 2) \leq 11$, but time out after twenty minutes when trying to determine if $ M(8, 3, 2) \in \{ 9, 10 \}$. Note that $ k=3, l=2$ is the least complex choice for the other parameters, so it seems like there’s not much hope to find $ M(n, k, l)$ for any seriously large parameters, like, say, $ k=6$.


These experiments suggest what SMT solvers can offer above and beyond ILP solvers. Disjunctions and implications are notoriously hard to model in an ILP. You often need to add additional special variables, or do tricks like the one I did that only work in some situations and which can mess with the efficiency of the solver. With SMT, implications are trivial to model, and natively supported by the solver.

Aside from reading everything I could find on Z3, there seems to be little advice on modeling to help the solver run faster. There is a ton of literature for this in ILP solvers, but if any readers see obvious problems with my SMT models, please chime in! I’d love to hear from you. Even without that, I am pretty impressed by how fast the solves finish for this subset cover problem (which this experiment has shown me is apparently a very hard problem).

However, there’s an elephant in the room. These models are all satisfiability/feasibility checks on a given solution. What is not tested here is optimization, in the sense of having the number of choice sets used be minimized directly by the solver. In a few experiments on even simpler models, z3 optimization is quite slow. And while I know how I’d model the ILP version of the optimization problem, given that it’s quite slow to find a feasible instance when the optimal number of sets is given as a parameter, it seems unlikely that it will be fast when asked to optimize. I will have to try that another time to be sure.

Also, I’d like to test the ILP models on Gurobi, but I don’t have a personal license. There’s also the possibility that I can come up with a much better ILP formulation, say, with a tighter LP relaxation. But these will have to wait for another time.

In the end, this experiment has given me some more food for thought, and concrete first-hand experience, on the use of SMT solvers.

Earthmover Distance

Problem: Compute distance between points with uncertain locations (given by samples, or differing observations, or clusters).

For example, if I have the following three “points” in the plane, as indicated by their colors, which is closer, blue to green, or blue to red?


It’s not obvious, and there are multiple factors at work: the red points have fewer samples, but we can be more certain about the position; the blue points are less certain, but the closest non-blue point to a blue point is green; and the green points are equally plausibly “close to red” and “close to blue.” The centers of masses of the three sample sets are close to an equilateral triangle. In our example the “points” don’t overlap, but of course they could. And in particular, there should probably be a nonzero distance between two points whose sample sets have the same center of mass, as below. The distance quantifies the uncertainty.


All this is to say that it’s not obvious how to define a distance measure that is consistent with perceptual ideas of what geometry and distance should be.

Solution (Earthmover distance): Treat each sample set $ A$ corresponding to a “point” as a discrete probability distribution, so that each sample $ x \in A$ has probability mass $ p_x = 1 / |A|$. The distance between $ A$ and $ B$ is the optional solution to the following linear program.

Each $ x \in A$ corresponds to a pile of dirt of height $ p_x$, and each $ y \in B$ corresponds to a hole of depth $ p_y$. The cost of moving a unit of dirt from $ x$ to $ y$ is the Euclidean distance $ d(x, y)$ between the points (or whatever hipster metric you want to use).

Let $ z_{x, y}$ be a real variable corresponding to an amount of dirt to move from $ x \in A$ to $ y \in B$, with cost $ d(x, y)$. Then the constraints are:

  • Each $ z_{x, y} \geq 0$, so dirt only moves from $ x$ to $ y$.
  • Every pile $ x \in A$ must vanish, i.e. for each fixed $ x \in A$, $ \sum_{y \in B} z_{x,y} = p_x$.
  • Likewise, every hole $ y \in B$ must be completely filled, i.e. $ \sum_{y \in B} z_{x,y} = p_y$.

The objective is to minimize the cost of doing this: $ \sum_{x, y \in A \times B} d(x, y) z_{x, y}$.

In python, using the ortools library (and leaving out a few docstrings and standard import statements, full code on Github):

from ortools.linear_solver import pywraplp

def earthmover_distance(p1, p2):
    dist1 = {x: count / len(p1) for (x, count) in Counter(p1).items()}
    dist2 = {x: count / len(p2) for (x, count) in Counter(p2).items()}
    solver = pywraplp.Solver('earthmover_distance', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

    variables = dict()

    # for each pile in dist1, the constraint that says all the dirt must leave this pile
    dirt_leaving_constraints = defaultdict(lambda: 0)

    # for each hole in dist2, the constraint that says this hole must be filled
    dirt_filling_constraints = defaultdict(lambda: 0)

    # the objective
    objective = solver.Objective()

    for (x, dirt_at_x) in dist1.items():
        for (y, capacity_of_y) in dist2.items():
            amount_to_move_x_y = solver.NumVar(0, solver.infinity(), 'z_{%s, %s}' % (x, y))
            variables[(x, y)] = amount_to_move_x_y
            dirt_leaving_constraints[x] += amount_to_move_x_y
            dirt_filling_constraints[y] += amount_to_move_x_y
            objective.SetCoefficient(amount_to_move_x_y, euclidean_distance(x, y))

    for x, linear_combination in dirt_leaving_constraints.items():
        solver.Add(linear_combination == dist1[x])

    for y, linear_combination in dirt_filling_constraints.items():
        solver.Add(linear_combination == dist2[y])

    status = solver.Solve()
    if status not in [solver.OPTIMAL, solver.FEASIBLE]:
        raise Exception('Unable to find feasible solution')

    return objective.Value()

Discussion: I’ve heard about this metric many times as a way to compare probability distributions. For example, it shows up in an influential paper about fairness in machine learning, and a few other CS theory papers related to distribution testing.

One might ask: why not use other measures of dissimilarity for probability distributions (Chi-squared statistic, Kullback-Leibler divergence, etc.)? One answer is that these other measures only give useful information for pairs of distributions with the same support. An example from a talk of Justin Solomon succinctly clarifies what Earthmover distance achieves

Screen Shot 2018-03-03 at 6.11.00 PM.png

Also, why not just model the samples using, say, a normal distribution, and then compute the distance based on the parameters of the distributions? That is possible, and in fact makes for a potentially more efficient technique, but you lose some information by doing this. Ignoring that your data might not be approximately normal (it might have some curvature), with Earthmover distance, you get point-by-point details about how each data point affects the outcome.

This kind of attention to detail can be very important in certain situations. One that I’ve been paying close attention to recently is the problem of studying gerrymandering from a mathematical perspective. Justin Solomon of MIT is a champion of the Earthmover distance (see his fascinating talk here for more, with slides) which is just one topic in a field called “optimal transport.”

This has the potential to be useful in redistricting because of the nature of the redistricting problem. As I wrote previously, discussions of redistricting are chock-full of geometry—or at least geometric-sounding language—and people are very concerned with the apparent “compactness” of a districting plan. But the underlying data used to perform redistricting isn’t very accurate. The people who build the maps don’t have precise data on voting habits, or even locations where people live. Census tracts might not be perfectly aligned, and data can just plain have errors and uncertainty in other respects. So the data that district-map-drawers care about is uncertain much like our point clouds. With a theory of geometry that accounts for uncertainty (and the Earthmover distance is the “distance” part of that), one can come up with more robust, better tools for redistricting.

Solomon’s website has a ton of resources about this, under the names of “optimal transport” and “Wasserstein metric,” and his work extends from computing distances to computing important geometric values like the barycenter, computational advantages like parallelism.

Others in the field have come up with transparency techniques to make it clearer how the Earthmover distance relates to the geometry of the underlying space. This one is particularly fun because the explanations result in a path traveled from the start to the finish, and by setting up the underlying metric in just such a way, you can watch the distribution navigate a maze to get to its target. I like to imagine tiny ants carrying all that dirt.

Screen Shot 2018-03-03 at 6.15.50 PM.png

Finally, work of Shirdhonkar and Jacobs provide approximation algorithms that allow linear-time computation, instead of the worst-case-cubic runtime of a linear solver.

Boolean Logic in Polynomials

Problem: Express a boolean logic formula using polynomials. I.e., if an input variable $ x$ is set to $ 0$, that is interpreted as false, while $ x=1$ is interpreted as true. The output of the polynomial should be 0 or 1 according to whether the formula is true or false as a whole.

Solution: You can do this using a single polynomial.

Illustrating with an example: the formula is $ \neg[(a \vee b) \wedge (\neg c \vee d)]$ also known as

not((a or b) and (not c or d))

The trick is to use multiplication for “and” and $ 1-x$ for “not.” So $ a \wedge b$ would be $ x_1 x_2$, and $ \neg z$ would be $ 1-z$. Indeed, if you have two binary variables $ x$ and $ y$ then $ xy$ is 1 precisely when both are 1, and zero when either variable is zero. Likewise, $ 1-x = 1$ if $ x$ is zero and zero if $ x$ is one.

Combine this with deMorgan’s rule to get any formula. $ a \vee b = \neg(\neg a \wedge \neg b)$ translates to $ 1 – (1-a)(1-b)$. For our example above,

$ \displaystyle f(x_1, x_2, x_3, x_4) = 1 – (1 – (1-a)(1-b))(1 – c(1-d))$

Which expands to

$ \displaystyle 1 – a – b + ab + (1-d)(ac + bc – abc)$

If you plug in $ a = 1, b = 0, c = 1, d = 0$ you get True in the original formula (because “not c or d” is False), and likewise the polynomial is

$ \displaystyle 1 – 1 – 0 + 0 + (1-0)(1 + 0 – 0) = 1$

You can verify the rest work yourself, using the following table as a guide:

0, 0, 0, 0 -&gt; 1
0, 0, 0, 1 -&gt; 1
0, 0, 1, 0 -&gt; 1
0, 0, 1, 1 -&gt; 1
0, 1, 0, 0 -&gt; 0
0, 1, 0, 1 -&gt; 0
0, 1, 1, 0 -&gt; 1
0, 1, 1, 1 -&gt; 0
1, 0, 0, 0 -&gt; 0
1, 0, 0, 1 -&gt; 0
1, 0, 1, 0 -&gt; 1
1, 0, 1, 1 -&gt; 0
1, 1, 0, 0 -&gt; 0
1, 1, 0, 1 -&gt; 0
1, 1, 1, 0 -&gt; 1
1, 1, 1, 1 -&gt; 0

Discussion: This trick is used all over CS theory to embed boolean logic within polynomials, and it makes the name “boolean algebra” obvious, because it’s just a subset of normal algebra.

Moreover, since boolean satisfiability—the problem of algorithmically determining if a boolean formula has a satisfying assignment (a choice of variables evaluating to true)—is NP-hard, this can be used to show certain problems relating to multivariable polynomials is also hard. For example, finding roots of multivariable polynomials (even if you knew nothing about algebraic geometry) is hard because you’d run into NP-hardness by simply considering the subset of polynomials coming from boolean formulas.

Here’s a more interesting example, related to the kinds of optimization problems that show up in modern machine learning. Say you want to optimize a polynomial $ f(x)$ subject to a set of quadratic equality constraints. This is NP-hard. Here’s why.

Let $ \varphi$ be a boolean formula, and $ f_\varphi$ its corresponding polynomial. First, each variable $ x_i$ used in the polynomial can be restricted to binary values via the constraint $ x_i(x_i – 1) = 0$.

You can even show NP-hardness if the target function to optimize is only quadratic. As an exercise, one can express the subset sum problem as a quadratic programming problem using similar choices for the constraints. According to this writeup you even express subset sum as a quadratic program with linear constraints.

The moral of the story is simply that multivariable polynomials can encode arbitrary boolean logic.