Checking consistency of a system of linear equations and inequalities

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.)

-
Its difficult to give an appropriate answer here but a very simple approach would be to use view this as a convex feasibility problem: You try to decide if the intersection of the solution space of the linear equation with the cube defined by the inequalities is not empty. Since your linear systems seem to be small and it looks like projecting a point onto the polytope formed by the linear constraints is also easy, alternating projections should be easy and fast. –  Dirk Mar 4 '12 at 10:43
Have a look at constraint-logic programming (in particular, over the reals -- CLP(R)) which is designed for just such problems. Typically these are implemented as Prolog libraries (Sicstus, SWI), but I believe that stand-alone versions are available too. –  J.J. Green Mar 4 '12 at 11:21
For these systems solving them is no harder than detecting feasibility. So just add an objective function, e.g., x+y+w+z, and now you have a linear program. Stuff this into an LP solver and it will tell if your feasible region is empty. If you read the manual, there may a flag you can set to avoid dding the objective function. This question would be just as well answered at math.stackexchange, and is not appropriate for this site. –  Chris Godsil Mar 4 '12 at 15:10
Chris, the inequalities there are strict, and these need a bit of work. A solution to your LP would likely violate the feasibility of these strict inequalities. –  Dima Pasechnik Mar 4 '12 at 15:28
Dima: you're right, of course. –  Chris Godsil Mar 4 '12 at 17:02

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 $$\label{eee} v\geq 0,\ M^\top v=0,\ v^\top b\leq 0.\qquad \qquad (*)$$

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.

-

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.

-