Find the Discrete Pair of {x,y} that Satisfy Inequality Constriants

Question

I have a few inequalities regarding {x,y}, that satisfies the following equations:

x>=0
y>=0
f(x,y)=x^2+y^2>=100
g(x,y)=x^2+y^2<=200

Note that x and y must be integer.

Graphically it can be represented as follows, the blue region is the region that satisfies the above inequalities:

The question now is, is there any function in Matlab that finds every admissible pair of {x,y}? If there is an algorithm to do this kind of thing I would be glad to hear about it as well.

Of course, one approach we can always use is brute force approach where we test every possible combination of {x,y} to see whether the inequalities are satisfied. But this is the last resort, because it's time consuming. I'm looking for a clever algorithm that does this, or in the best case, an existing library that I can use straight-away.

The x^2+y^2>=100 and x^2+y^2<=200 are just examples; in reality f and g can be any polynomial functions of any degree.

Edit: C# code are welcomed as well.

Dom

Answer 1

This is surely not possible to do in general for a general set of polynomial inequalities, by any method other than enumerative search, even if there are a finite number of solutions. (Perhaps I should say not trivial, as it is possible. Enumerative search will work, subject to floating point issues.) Note that the domain of interest need not be simply connected for higher order inequalities.

Edit: The OP has asked about how one might proceed to do a search.

Consider the problem

x^3 + y^3 >= 1e12
x^4 + y^4 <= 1e16

x >= 0, y >= 0

Solve for all integer solutions of this system. Note that integer programming in ANY form will not suffice here, since ALL integer solutions are requested.

Use of meshgrid here would force us to look at points in the domain (0:10000)X(0:10000). So it would force us to sample a set of 1e8 points, testing every point to see if they satisfy the constraints.

A simple loop can potentially be more efficient than that, although it will still require some effort.

% Note that I will store these in a cell array,
% since I cannot preallocate the results.
tic
xmax = 10000;
xy = cell(1,xmax);
for x = 0:xmax
  % solve for y, given x. This requires us to
  % solve for those values of y such that
  %   y^3 >= 1e12 - x.^3
  %   y^4 <= 1e16 - x.^4
  % These are simple expressions to solve for.
  y = ceil((1e12 - x.^3).^(1/3)):floor((1e16 - x.^4).^0.25);
  n = numel(y);
  if n > 0
    xy{x+1} = [repmat(x,1,n);y];
  end
end
% flatten the cell array
xy = cell2mat(xy);
toc

The time required was...

Elapsed time is 0.600419 seconds.

Of the 100020001 combinations that we might have tested for, how many solutions did we find?

size(xy)
ans =
           2     4371264

Admittedly, the exhaustive search is simpler to write.

tic
[x,y] = meshgrid(0:10000);
k = (x.^3 + y.^3 >= 1e12) & (x.^4 + y.^4 <= 1e16);
xy = [x(k),y(k)];
toc

I ran this on a 64 bit machine, with 8 gig of ram. But even so the test itself was a CPU hog.

Elapsed time is 50.182385 seconds.

Note that floating point considerations will sometimes cause a different number of points to be found, depending on how the computations are done.

Finally, if your constraint equations are more complex, you might need to use roots in the expression for the bounds on y, to help identify where the constraints are satisfied. The nice thing here is it still works for more complicated polynomial bounds.

Mut

@woodchips, it's not possible, and there is no algorithm as well? – Dom

@Ngu: I think that one of the implications of a problem being impossible is that there cannot be an algorithm to solve it. I agree with @woodchips. If you doubt this, start writing down, with pen and paper, the first few (x,y) pairs that satisfy your inequalities. Then stop for a while and think. – Aiken

For the general problem, NO, there is no simple solution. You can probably do no better than to set up a loop on one variable, say x. Fix x at some integer value, then solve for all y that satisfy the constraints. A problem is, it may take some effort just show that once you get to a point that there are no solutions for y for a given x, that there is not a larger value of x that might yet yield further solutions. Interval arithmetic might help on this aspect. – Mut

@woodchips, what if the polynomial is known? – Dom

@High, what if the polynomial is known? Is it still unsolvable? – Dom

Simply knowing the polynomial coefficients does not make it easy. Even in the case of linear constraints, enumeration of all integer solutions is not trivial. I am NOT saying it is not solvable, but that you are stuck with a loop. Again, brute force can be viewed as a loop on x, then solving for all possible integer solutions in y, given that value of x. The biggest problem here is knowing when to stop with certainty. That will be difficult for high order polynomials. – Mut

In the case of linear constraints, this is Linear Integer Programming, which is NP-complete but has practical solutions for small instances, such as the cutting-plane method. – Kannada

@Jouni: No. Read the question as asked. The OP asked to find ALL solutions. LP will yield A solution. LP does not find the entire set of solutions. – Mut

@woodchips, sure, I should have said "linear constraints and even deciding if there are any solutions". There can be infinitely many solutions, and even for bounded programs an exponential number of corner point solutions. Still, in small cases you may be able to enumerate them. – Kannada

@woodchips, You can probably do no better than to set up a loop on one variable, say x.... not sure what you mean here, are you saying that setup a loop on one variable and then solve it, is no better than brute enumerative search? Why this is so? – Dom

y = ceil((1e12 - x.^3).^(1/3)):floor((1e16 - x.^4).^0.25); I think that's a clever one. All the in the range y were collected at once, instead of the need to compute the y value one by one and rejecting the ones that dont fit. – Dom

Admittedly, this was a pretty solution. It works nicely for a pair of constraints that have simple solutions once x is given. However, if you have more general 2 variable polynomial expressions, use roots on the result after substituting x. With a little play, that will yield intervals in y that must work for a given x. – Mut