# Find a solution of large system of inequalities

I have a large system of homogenous inequalities involving 33 real unknowns of the form

$$\vec{F}(z_i)^T \cdot \vec{X}>0\,$$

where $$\vec{X} = \left(x_1,...,x_{24}\right)^T$$ are the unknowns and $$\vec{F}(z_i)$$ are functions depending on a real parameter $$|z_i|<1$$ which parametrises the system, that is to each value of $$z_i$$ is associated an inequality. An example might be $$\vec{F}(z) = \left(z^2+z\,,\, z^3+1/z,...,0.1 z^4 + 1/\sqrt{z}\right)^T$$ but I am interested into more general cases.

I am interested in finding a point $$\vec{X}_0$$ which satisfy the system of inequalities resulting from the very big set of parameters $$\{z_i\}$$. Let say, we select 100 different parameters $$z_i$$, then we have 100 inequalities to solve.

I have tried to solve this problem with Mathematica, which has a function FindInstance which does this job. This function works very well for smaller systems, whereas it fails for large system of inequalities.

Questions

• Is there a way to study this problem by using Python or C++?

• Are there already dedicated libraries to solve such kind of problems?

• I want to find a set of points $$\{\vec{X}\}$$ which solves the system.

• Your problem seems to be finding a feasible point for a set of inequality constraints. Mar 7 '19 at 13:36
• @nicoguaro actually it is a feasible test problem. I am not an expert of this, and I am reading something about linear programming around in the internet. Any suggestion, is welcome! Mar 7 '19 at 13:37
• Have you checked SciPy.optimize? Mar 7 '19 at 15:22
• Are you really interested in the intersection of open half spaces? That's numerically not useful. Mar 9 '19 at 1:21
• @WolfgangBangerth yes, I am really interested in the strictly positive inequality. Do you mean that's hard to solve this problem numerically? I understand the difficulties... Mar 9 '19 at 14:45

The goal is to compute an $$x\in\mathbb{R}^{n}$$ to satisfy the strict inequalities $$f(z_{i})^{T}x>0\qquad\forall z_{i}\in\{1,\ldots,m\}.\tag{1}$$ If we write $$\epsilon>0$$ as the smallest margin of feasibility, that is $$\epsilon=\min_{i}\{f(z_{i})^{T}x\},$$ then (1) is equivalent to $$f(z_{i})^{T}x\ge\epsilon>0\qquad\forall z_{i}\in\{1,\ldots,m\}.\tag{2}$$ Proposition. There exists a pair $$(x,\epsilon)$$ satisfying the strict inequality (2) if and only if there exists $$y\in\mathbb{R}^{n}$$ satisfying the nonstrict inequality (3) below $$f(z_{i})^{T}y\ge1\qquad\forall z_{i}\in\{1,\ldots,m\}.\tag{3}$$

Proof. If. Any $$y$$ satisfying (3) also yields the pair $$(y,1)$$ satisfying (2). Only if. By contradiction, suppose that there exists no $$y$$ satisfying (3) but there exists a pair $$(x,\epsilon)$$ satsifying (2). Then $$y=x/\epsilon$$ satisfies (3)--- a contradiction.

Hence, we can just go and find a choice of $$y$$ satisfying (3). Any standard off-the-shelf linear programming solver will solve (3) in standard canonical form $$\begin{gather*} \min_{z}\quad c^{T}z\quad\text{subject to }\quad Az=b,\quad z\ge0,\\ \max_{y}\quad b^{T}y\quad\text{subject to }\quad A^{T}y\le c, \end{gather*}$$ with data $$A=-\begin{bmatrix}f(z_{1}) & \cdots & f(z_{m})\end{bmatrix}^T,\qquad b=0,\qquad c=-\begin{bmatrix}1 & \cdots & 1\end{bmatrix}^T.$$ A problem with $$n\approx m\approx100$$ is extremely small. In this case, you can also try 3 lines of MATLAB code using ADMM. (See the Boyd et al. survey for more details. At the risk of a little self-promotion, the specific derivations can be found in Section IV-C of our tutorial paper): \begin{align*} y^{k+1} & =(AA^{T})^{-1}A\left[c-s^{k}-z^{k}/t\right]\\ s^{k+1} & =\max\{0,c-A^{T}y^{k+1}-z^{k}/t\}\\ z^{k+1} & =z^{k}+t(A^{T}y^{k+1}+s^{k+1}-c) \end{align*} Here, you can pick any step-size $$t>0$$ and just set the initial points to zero. See also Jacob Mattingley's page on ridiculously short LP solvers.

• Thank you for this answer. I will try to learn how to use linear programming solvers and will come back to you Mar 13 '19 at 8:31
• I should note that if a choice of $x$ does not exist, then you can use Lagrange duality to certify it. The geometrical interpretation is to generate a separating hyperplane that has the space $\mathrm{span}(A)$ on one side and the cone $x>0$ on the other. Such a plane is indeed naturally computed by any LP solver Mar 13 '19 at 16:15

Here are two considerations that might help you figure out whether a solution exists:

• First, if there is a vector $$X$$ so that $$F(z_i) \cdot X > 0$$, then all vectors $$Y=\alpha X$$ for any choice of $$\alpha>0$$ are also solutions because $$F(z_i) \cdot Y = \alpha \underbrace{(F(z_i)\cdot X)}_{\ge 0} > 0$$. In other words, if a solution exists, then there are in fact infinitely many solutions.

• Second, let's forget that the $$F(z_i)$$ vectors are really some kind of parameterized version. They are really just a bunch of vectors $$F_i$$ that happen to be created in a special way. Then the inequality $$F_i \cdot X>0$$ really just means that $$X$$ must lie in the (open) half space cut out of all of $${\mathbb R}^n$$ by the plane that is perpendicular to the vector $$F_i$$. A vector $$X$$ that satisfies all of the inequalities $$F_i \cdot X>0, i=1,\ldots,n$$ must therefore lie in the intersection of the $$n$$ half spaces. The question of existence of a solution then boils down to the question whether the intersection of the $$n$$ half spaces is non-empty. This is a geometric question. For example, if you were in 2d, then if you had vectors $$F_1=(1,0), F_2=(-1,0)$$ then the (open) half-spaces are mutually exclusive and so no solution exists. The same would be true if you had $$F_1=(1,0), F_2=(-1,1), F_3=(-1,-1)$$, even though the intersection of each pair of half spaces is non-empty.

Out of this second consideration, I suspect that you can build an algorithm that determines whether the intersection of the first $$k$$ half spaces with the $$(k+1)$$st half space is non-empty. I'd have to think about this in more depth to come up with an algorithm, but maybe you can do a literature search about the intersection of half spaces to come up with something useful. The point worth mentioning is that the intersection of your half spaces is a cone rooted at the origin that I suspect can be described by the intersection of the dividing planes of each half space; these are low dimensional objects easily amenable to to linear algebra.

Of course, there is also the possibility to use the special structure of the $$F_i$$ to prove that a solution does or does not exist. For example, if you can show that there is a direction $$Y$$ so that all of the $$F_i=F(z_i)$$ have an angle less than 90 degree from $$Y$$ (i.e., $$Y\cdot F_i>0$$), then it is geometrically clear that $$Y$$ is a vector that satisfies the desired inequalities. Of course, you may not know $$Y$$, but maybe the statement is true for $$Y=F(z_1)$$, for example.

• Thank you for this answer. I have already thought the geometrical interpretation of the problem, which of course coincides with yours. For instance, I have tried to construct the convex-hull generated by the different points $\vec{F}_i$ and see whether or not $\vec{0}$ belongs to the convex-hull. This might be a proof that the system in impossible. However, constructing high-dimensional convex hulls is challenging and it seems to me that $\vec{0}$ belongs likely to the boundary of the convex-hull. So, since $0$ is the same as $10^{-30}$, this method might not be appropriate. Mar 13 '19 at 8:37
• I think the convex hull is not a useful tool here. That's because you can scale these $F_i$ arbitrarily: If you replace once constraint $F_i \cdot X>0$ by $(2F_i)\cdot X>0$, nothing actually changes. So you need to find a representation of your problem that takes this into account. The representation as a cone I mentioned is something that respects this invariance. Mar 13 '19 at 21:05

Since $$\vec{F}(z_i)$$ is a constant once the parameter $$z_i$$ has been chosen, the problem amounts to solving a system $$\mathbf{A}\vec x>0$$ for which $$\vec x =0$$ is trivially a feasible point.

EDIT

The problem is that no optimization software implements strict inequalities. One reason for this is that the difference between say, 1e-999 and 0 is really, really small. So, operationally, a strict equality and an equality shakeout to the same answer. Let's explore this a bit.

cvxpy is a Python package that allows you to write a linear program in a generic way and then pass it to one of a large number solvers. Like most such programs, it does not allow for strict inequalities.

I've written a representation of your problem below:

#!/usr/bin/env python3
import cvxpy as cp
import numpy as np

def f(z):
return np.array([
z**2+z-3,
z**3+1/z,
0.1*z**4+1/np.sqrt(z),
-4
])

#Generate 100 random numbers in [-1,1)
zs = np.random.uniform(low=-1, high=1, size=100)
#Run the function on each number
zs = [f(z) for z in zs]
#Strip out results with a NaN
zs = [z for z in zs if not any(np.isnan(z))]

zs = [[-1,-1,-1,-1],[1,1,1,1]]

#x vector of length four to match length of vector from f
x = cp.Variable(4)

#Multiply f(z)'s by x's and build constraints
cons = [z*x>=0 for z in zs]

#Constant objective implies a problem in which we only want to find a feasible
#point
obj = cp.Maximize(1)

#Create a problem with the objective and constraints
prob = cp.Problem(obj, cons)

#Solve problem, get optimal value
val = prob.solve()

if val==-np.inf:
print("NO SOLUTION FOUND")
else:
print(x.value)


If we run the above code, we find that the solver indeed finds the 0-vector as an answer.

Now, you would like your answer to be strictly greater than zero. Since there is no solver that will do this (that I know of), one way to achieve it is to ask the solver to find a solution to a different problem. Rather than solving $$Ax>0$$, we can ask the solver to find $$Ax\ge\epsilon$$. Where $$\epsilon$$ is some small value. Since all you want is a feasible point, this is a reasonable approach since any solution to $$Ax\ge\epsilon$$ is also a solution to $$Ax>0$$. This is the simple solution you thought you were looking for.

We can implement this by changing the line:

cons = [z*x>=0 for z in zs]


to

cons = [z*x>=0.0001 for z in zs]


Running the code again, I get the answer:

[ 1.34351053e-09 -9.62045447e-10  5.11738862e-09 -3.99990907e-05]


(You may get a different answer because the solvers leverage stochasticity.) This looks promising! But there's a caveat here... What if we use this program?

#!/usr/bin/env python3
import cvxpy as cp
import numpy as np

#A system with no solution
zs = [[-1,-1,-1,-1],[1,1,1,1]]

#x vector of length four to match length of vector from f
x = cp.Variable(4)

#Multiply f(z)'s by x's and build constraints
cons = [z*x>=0.0001 for z in zs]

#Constant objective implies a problem in which we only want to find a feasible
#point
obj = cp.Maximize(1)

#Create a problem with the objective and constraints
prob = cp.Problem(obj, cons)

#Solve problem, get optimal value
val = prob.solve()

if val==-np.inf:
print("NO SOLUTION FOUND")
else:
print(x.value)


It's obvious from inspection that there can be no solution satisfying this system. Interestingly, however, when you run the program (using the default solver) the 0-vector is returned as an answer! If we modify the program to read:

cons = [z*x>=0.001 for z in zs]


the same thing happens. Only when we get to:

cons = [z*x>=0.01 for z in zs]


do we finally get the correct response: that there is no solution.

There are a few reasons for this:

• Internally, the solver is using a floating-point representation whose limited precision results in it thinking it's solved the problem when it really hasn't. (You could deal with this by using a rational-number solver that does its calculations using fractional representations.)
• More generally, the solver may not be numerically robust. The subfield of "robust optimization" can be leveraged to find disciplined ways of handling this.
• Philosophically, asking the solver to differentiate between small values of epsilon and 0 is silly. Say you're optimizing the floor plan of a house using numbers which represent metres and choose epsilon as 1e-10. You're asking for a solution that differs from zero by the width of an atom. Say you're calculating a solar trajectory with numbers representing astronomical units (1 AU is the distance from Earth to the sun - 93 million miles): the difference between 1e-10 and 0 is 50 feet (the width of a house).

Perhaps the simplest way of dealing with the problems above is to rescale your system so that small values are unimportant. For instance, rather than measuring my floorplan in meters, I could measure it in millimeters.

It's worth noting that there is an entire class of problems for which your question has trivial and reliable solutions. If any column of your $$A$$ matrix contains only positive values greater than zero, then setting that column's corresponding $$x$$ value to 1 and all other $$x$$ values to 0 provides an answer. Similarly, if any column contains only negative values then choosing -1 for the corresponding $$x$$ value and 0 elsewhere provides an answer.

• Why $\vec{x}=0$? Mar 8 '19 at 6:25
• Because if you multiply any $\vec F(z_i)$ by $\vec 0$, you get 0, which satisfies your constraint (since you use a strict inequality, in practice you perturb $\vec 0$ by an infinitesimal positive amount). Mar 8 '19 at 7:42
• I don't get this answer. The inequality is strict: $Ax>0$. $x=0$ does not satisfy the inequality. Mar 9 '19 at 4:58
• @WolfgangBangerth: Edited. Mar 12 '19 at 0:56
• It still doesn't make sense. For example, assume $x\in\mathbb R$ and that $A$ is just the matrix with entries $[1 ; -1]$, corresponding to the inequalities $x>0$ and $x<0$. The problem has no feasible point. It also does not have a feasible point if you replace the zero by epsilon: $x>\varepsilon$ and $-x>\varepsilon$ still has no solution. But the problem with $\ge 0$ does have a feasible point, namely $x=0$. Mar 12 '19 at 3:57