2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Your problem seems to be finding a feasible point for a set of inequality constraints. $\endgroup$ – nicoguaro Mar 7 at 13:36
  • $\begingroup$ @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! $\endgroup$ – apt45 Mar 7 at 13:37
  • $\begingroup$ Have you checked SciPy.optimize? $\endgroup$ – nicoguaro Mar 7 at 15:22
  • $\begingroup$ Are you really interested in the intersection of open half spaces? That's numerically not useful. $\endgroup$ – Wolfgang Bangerth Mar 9 at 1:21
  • $\begingroup$ @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... $\endgroup$ – apt45 Mar 9 at 14:45
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for this answer. I will try to learn how to use linear programming solvers and will come back to you $\endgroup$ – apt45 Mar 13 at 8:31
  • $\begingroup$ 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 $\endgroup$ – Richard Zhang Mar 13 at 16:15
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – apt45 Mar 13 at 8:37
  • $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Mar 13 at 21:05
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Why $\vec{x}=0$? $\endgroup$ – apt45 Mar 8 at 6:25
  • $\begingroup$ 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). $\endgroup$ – Richard Mar 8 at 7:42
  • 3
    $\begingroup$ I don't get this answer. The inequality is strict: $Ax>0$. $x=0$ does not satisfy the inequality. $\endgroup$ – Wolfgang Bangerth Mar 9 at 4:58
  • $\begingroup$ @WolfgangBangerth: Edited. $\endgroup$ – Richard Mar 12 at 0:56
  • 1
    $\begingroup$ 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$. $\endgroup$ – Wolfgang Bangerth Mar 12 at 3:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.