Counting solutions to linear inequalities [closed]

Question

I am trying to solve a problem which I have reduced down to counting the number of integer solutions to a number of linear inequalities. I need to be able to count the number of solutions for any number of variables c_1, ..., c_n, but for n=3 the equations could be written as:

The equations. http://silicon.appspot.com/readdoc?id=155604

Now, I know the values of n and r in advance and wish to find the number of (c_1, ..., c_n) solutions that exist.

Can this be done efficiently (faster than enumerating the solutions)? (If so: how?; if not: why?)

Answer 1

Suppose you had some code to produce all the solutions.

(For the parameter z here, pass 9. It's the number on the right-hand side of the inequalities. Note that this code only works when r is positive.)

from math import floor, ceil

def iter_solutions(r, n, z):
    c = [None] * n
    def iter_solutions_bounded(k, pick):
        # pick is the last pick, if any, and 0 otherwise
        assert (1 <= k < n and pick == c[k]) or (k == n and pick == 0)

        min_ck = int(ceil(-pick / r))
        max_ck = int(floor((z - pick) / r))
        if k == 1:
            for ck in range(max(min_ck, 0), min(max_ck, z) + 1):
                c[0] = ck
                yield c
        else:
            for ck in range(min_ck, max_ck + 1):
                c[k - 1] = ck
                for soln in iter_solutions_bounded(k - 1, ck):
                    yield soln
    return iter_solutions_bounded(n, 0)

You can convert this into code that merely counts the solutions simply by deleting all the code that refers to c and adding up the number of solutions that would have been yielded. Finally, you can improve the performance by memoization.

from math import floor, ceil

def memoize(f):
    cache = {}
    def g(*args):
        if args in cache:
            return cache[args]
        tmp = cache[args] = f(*args)
        return tmp
    return g

def len_range(a, b):
    if a <= b:
        return b - a
    return 0

def count_solutions(r, n, z):
    @memoize
    def count_solutions_bounded(k, pick):
        min_ck = int(ceil(-pick / r))
        max_ck = int(floor((z - pick) / r))
        if k == 1:
            return len_range(max(min_ck, 0), min(max_ck, z) + 1)
        else:
            return sum(count_solutions_bounded(k - 1, ck) for ck in range(min_ck, max_ck + 1))
    return count_solutions_bounded(n, 0)

Some possible improvements:

• If it's true that c₁ ... c_n are always ≤ z, then detecting that and immediately returning 0 would help a lot for large n. In fact it would reduce the running time to a lightning-fast O(nz).

• If it is intended that c₁ ... c_n all be non-negative, that's even better. Making the appropriate changes to min_ck and max_ck would make this O(nz) with a smaller constant, and the cache could be a flat 2D array instead of the slower hash table I've got.

• You might be able to do better by building the cache systematically, rather than populating it "on demand" the way this memoization code does. First build the whole cache for n=1, then for n=2, and so on. That way you could avoid recursion, and at each step you could throw away the cached data you no longer need (after calculating the results for n=2, you don't need the entries for n=1 anymore).

Answer 2

To solve this problem I would probably go into the realms of constraint programming. It seems like you have a classic all different constraint (a bit like the N-Queens problem). Have a go at one of the free constraint solvers listed below. That will give you a solution quite efficiently. It'll basically generate the whole search tree but with the nice All-Different constraint implementations out there, the tree will end up being pruned almost to nothing.

http://www.gecode.org/ http://minion.sourceforge.net/ http://jacop.osolpro.com/ http://research.microsoft.com/apps/pubs/default.aspx?id=64335

Here's the wikipedia list:
http://en.wikipedia.org/wiki/Constraint_programming#Constraint_programming_libraries_for_imperative_programming_languages

Answer 3

This is not really a full solution to your problem, but I think it may help or at least give you some ideas.

Your requirement that the solutions be integer makes this an NP problem. If we first consider the relaxation of the problem so that the domain is the real numbers, you are asking to solve the satisfiability problem of 0 <= A*c <= 1, where A is a matrix and c is your vector of unknowns. This is a standard linear program (LP with trivial objective), and can be solved efficiently (in polynomial time). You may want to use this as a first pass test to determine feasibility since if the relaxed LP has no solutions, your integer LP certainly has no solutions. A good LP solver will also return a feasible point if possible, and you may be able to round the vector's entries to find an integer solution.

Answer 4

As others have mentioned, if you wanted to maximize a linear objective function based on these constraints then you would have an integer linear programming problem, for which no efficient general solution exists. Instead you seem to be asking for the number of points in the feasible region, which is a different problem, but it too is complicated by having to have integer solutions.

The best idea I can come up with involves finding the points on the boundary of the feasible region, and using those to determine the number of points inside. This works well at accelerating "count the lattice points" type problems in lower dimensions, but the boundary is still just one dimension smaller than the volume in question. If your problem gets over a few dimensions then the problem will still be intractable, even if it is faster than enumerating all the solutions.