Question: 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, 0 0, 1, 1, 1  So 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.

Answer: 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  rather than using sparse modelling, you just loop don't use LP at all fix each variable at each of its extrema, and iteratively column-eliminate from the left-hand side attempt a least-squares solution of the linear system, and infer a high residual to mean that there is no solution  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.]