Question

I have a lot of systems of equations and inequalities of the following form:

$$ a_{1,1}x+a_{1,2}y+a_{1,3}z+a_{1,4}w = 2 $$ $$ \ldots $$ $$ 0 < x < 2 $$ $$ 0 < y < 2 $$ $$ 0 < z < 2 $$ $$ 0 < w < 2 $$

There are always at least two equations, and I probably won't consider cases with more than twenty equations. All coefficients $a_{i,j}$ are positive integers and some can be zero. We also have the property that $\sum_{j=1}^4a_{i,j}\geq3$ for all $i$. The solutions are real numbers.

I don't need to solve these systems, but I need to be able to tell whether there exists a solutions. If it isn't possible to tell for each system whether it is consistent or not, any method which identifies as many inconsistent systems as possible is greatly appreciated.

I have a few hundred millions of these systems, so I'm specifically looking for things that can easily be turned into a program. (I know the basic techniques to do this by hand, and am looking for some handy tricks that can be done by a computer. I have some programming experience, but not really with programming this kind of problems.)

Answer 1

There is a criterion for solvability of a system of strict inequalities $Mt\lt b$ due to Carver (cf. A.Schrijver "Theory of linear and integer programming", Sect. 3.7.8). It says that $Mt\lt b$ is solvable if and only if $v=0$ is the only solution of the system \begin{equation}\label{eee} v\geq 0,\ M^\top v=0,\ v^\top b\leq 0.\qquad \qquad (*)\end{equation}

Let us see how to get $Mt\lt b$. To do this, let $\zeta=\frac{1}{2}(x,y,z,w)$, and write your linear equations as $A\zeta=e$, where $e$ denotes all-1 vector. If this system has no solution, done. If it has only one solution, you can check it with your inequalities directly. If there are several solutions, linear algebra software will be able to rewrite your system in the form $(I\ B)\zeta'=d$, where $I$ is the identity matrix of size 1, 2, or 3, $B$ is a matrix of the appropriate size, and $\zeta'$ is a permutation of the original variables $\zeta$. In other words it gives you expressions $\zeta'_k=d_k-\sum_{j\neq k} B_{kj}\zeta'_j$, for $1\leq k\leq m$, and $m$ being 1, 2, or 3, depending upon $A$.

This reduces your original system to the system of strict inequalities in the remaining unexpressed $\zeta'$. (i.e. in $4-m$ variables).

Finally, apply Carver's criterion by solving a linear programming problem: $\max\sum_{i} v_i$ subject to $(*)$. If this maximum is strictly bigger than 0 then the original system $Mt\lt b$ has no solution, otherwise it does have one.

Answer 2

To simplify the notation, let $A$ be the coefficient matrix in a given instance of your problem, let $\xi = ((x,y,z,w)^T)/2$, and let $O = (0,0,0,0)^T$, $e = (1,1,1,1,...)^T$, $e_4 = (1,1,1,1)^T$. The problem then can be written in shorthand as $$ A \xi = e, O < \xi < e_4 $$ where $O < \xi < e_4$ is understood component wise.

To check if there is a feasible solution, proceed in two steps:

1. Check if there is a feasible solution of the linear system $A \xi = e$. If there is one, it can be found as $\xi = (A^TA)A^Te$ and therefore $A(A^TA)^{-1}A^Te = e$ must hold. In practice, compute the $QR$ decomposition of $A$, $A = QR$ where $R$ is square and upper triangular and $Q$ has the same dimensions as $A$ and satisfies $Q^TQ = I$ (identity matrix) and look at the system $R \xi = Q^T e$. If $R$ has full rank, you can find $\xi$ and compare $A \xi $ to $e$. If $R$ does not have full rank (i.e. there are zero rows at the bottom), this also tells you if there is a solution.

2. Suppose now there is a nontrivial solution of $A \xi = e$. Then choose a small number $\epsilon$, e.g. $\epsilon = 10^{-8}$, and solve the linear program $$ A \xi = e, \epsilon e_4 \le \xi \le (1 - \epsilon) e_4, c^T \xi \to \max $$ with any vector $c$. If there is a feasible solution to the full problem, it will show up as the optimum (and some of its components may be equal to $\epsilon$ or $1 - \epsilon$). By varying $\epsilon$ for these problems, you may be able to find solutions which are further in the interior of the four-dimensional cube in which your solution is supposed to be.

The entire method should be easily implementable in e.g. R (package lpSolve) and it is obvious how to parallelize it.