I have a system of equations where each equation is a linear equation with boolean constraints. For example:
x1 + x2 + x3 = 2
x1 + x4 = 1
x2 + x1 = 1
And each x_i is either 0 or 1. Sometimes there might be a small positive (<5) coefficient (for example x1 + 2 * x3 + x4 = 3. Basically a standard linear programming task. What I need to do is to find all x_i which are guaranteed to be 0 and all x_j which are guaranteed to be 1. Sorry if my terminology is not correct here but by guaranteed I mean that if you generate all possible solutions you in all of them all x_i will be 0 and in all of them x_j will be 1.
For example my equation has only 2 solutions:
1, 0, 1, 00, 1, 1, 1So you do not have guaranteed 0 and have x_3 as a guaranteed 1.
I know how to solve this problem with or-tools by generating all solutions and it works for my usecases (equations are pretty constrained so usually there are < 500 solutions although the number of variables is big enough to make the whole combinatorial search impossible).
The big problem is that I can't use that library (system restrictions above my control) and only libraries available in my case are numpy and scipy. I found that scipy has scipy.optimize.linprog.
It seems like I have found a way to generate one solution
import numpy as np
from scipy.optimize import linprog
A_eq = np.array([
[1, 1, 1, 0], # x1 + x2 + x3 = 2
[1, 0, 0, 1], # x1 + x4 = 1
[1, 1, 0, 0] # x1 + x2 = 1
])
b_eq = np.array([2, 1, 1])
c = np.zeros(4)
bounds = [(0, 1)] * 4
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs-ipm')
if res.success:
print(res.x)
But I can't find a way to generate all solutions. Also I am not sure whether there is a better way to do it as all I need to know is to find guaranteed values
P.S. this problem is important to me. I guarantee to add a 500 bounty on it, but system prevents me from doing it until 2 days will pass.
You don't need to (fully) brute-force, and you don't need to find all of your solutions. You just need to find solutions for which each of your variables meets each of their extrema. The following is a fairly brain-off LP approach with 2n² columns and 2mn rows. It's sparse, and for your inputs does not need to be integral. That said, I somewhat doubt it will be the most efficient method possible.
import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
import scipy.sparse as sp
lhs = np.array((
(1, 1, 1, 0),
(1, 0, 0, 1),
(1, 1, 0, 0),
))
rhs = np.array((2, 1, 1))
m, n = lhs.shape
# Variables: n * 2 (minimize, maximize) * n
c = sp.kron(
sp.eye_array(n),
np.array((
(+1,),
(-1,),
)),
)
b = np.tile(rhs, 2*n)
system_constraint = LinearConstraint(
A=sp.kron(sp.eye_array(2*n), lhs, format='csc'),
lb=b, ub=b,
)
result = milp(
c=c.toarray().ravel(), # must be dense
integrality=0,
bounds=Bounds(lb=0, ub=1),
constraints=system_constraint,
)
assert result.success
extrema = result.x.reshape((n, 2, n))
mins = extrema[:, 0]
maxs = extrema[:, 1]
vmins = np.diag(mins)
vmaxs = np.diag(maxs)
print('Solutions for minima on the diagonal:')
print(mins)
print('Solutions for maxima on the diagonal:')
print(maxs)
print('Variable minima:', vmins)
print('Variable maxima:', vmaxs)
print('Guaranteed 0:', vmaxs < 0.5)
print('Guaranteed 1:', vmins > 0.5)
Solutions for minima on the diagonal:
[[-0. 1. 1. 1.]
[ 1. 0. 1. -0.]
[ 1. 0. 1. -0.]
[ 1. 0. 1. -0.]]
Solutions for maxima on the diagonal:
[[ 1. 0. 1. -0.]
[-0. 1. 1. 1.]
[ 1. 0. 1. -0.]
[-0. 1. 1. 1.]]
Variable minima: [-0. 0. 1. -0.]
Variable maxima: [1. 1. 1. 1.]
Guaranteed 0: [False False False False]
Guaranteed 1: [False False True False]
There is a variant on this idea where
This somewhat naively assumes that all solutions will see integer values, and (unlike milp) does not have the option to set integrality=1. For demonstration I was forced to add a row to get a residual.
import numpy as np
lhs = np.array((
(1, 1, 1, 0),
(1, 0, 0, 1),
(1, 1, 0, 0),
(0, 0, 1, 1),
))
rhs = np.array((2, 1, 1, 1))
m, n = lhs.shape
epsilon = 1e-12
lhs_select = np.ones(n, dtype=bool)
for i in range(n):
lhs_select[i] = False
x0, (residual,), rank, singular = np.linalg.lstsq(lhs[:, lhs_select], rhs)
zero_solves = residual < epsilon
x1, (residual,), rank, singular = np.linalg.lstsq(lhs[:, lhs_select], rhs - lhs[:, i])
one_solves = residual < epsilon
lhs_select[i] = True
if zero_solves and not one_solves:
print(f'x{i}=0, solution {x0.round(12)}')
elif one_solves and not zero_solves:
print(f'x{i}=1, solution {x1.round(12)}')
x0=1, solution [-0. 1. 0.]
x1=0, solution [ 1. 1. -0.]
x2=1, solution [1. 0. 0.]
x3=0, solution [ 1. -0. 1.]