I'm trying to establish an upper bound on a fairly elaborate equation that comes from a tree. It has a small number of of inputs that can be encoded with "one-hot" variables. All numbers are integers. All terms are sums or products of other terms and booleans. Check out this question for details.
Following a suggestion in the response to that question, I've formulated my problem using Integer Linear Programming (ILP). I've encoded the problem with Google's OR-Tools in an SMT file to feed to Microsoft's Z3. For comparison, I also implemented a brute force solution in Python.
Much to my surprise, the brute force solution in Python was by far the fastest! I would have assumed it would be by far the slowest, so what am I missing?
Here's what the Python looks like (full code):
def w(c0, c1, c2, c3, c4, c5):
a0,e0,i0,o0,u0 = c0 # exactly one of these will be True
c1,h1,k1,m1,p1,t1 = c1 # (same for all the cN variables)
a2,e2,i2,o2,u2 = c2 # ...
l3,n3,r3,s3,y3 = c3 # ...
b4,d4,f4,g4,j4,v4,w4,x4,z4 = c4 # ...
c5,h5,k5,m5,p5,t5 = c5
n_0=1+d4
n_1=(s3*2)+(y3*2)
n_2=1+n_1
n_3=(d4*n_2)
n_4=(e2*n_0)+(i2*n_3)
n_5=(d4*2)
# ...
n_3293=n_3220+n_3281+n_3292
n_3294=(c5*n_2914)+(h5*n_2973)+(k5*n_3016)+(m5*n_3119)+(p5*n_3202)+(t5*n_3293)
n_3295=n_384+n_1082+n_1412+n_1938+n_2788+n_3294
return n_3295
Here's what the OR-Tools formulation looks like (full code):
from ortools.sat.python import cp_model
model = cp_model.CpModel()
cells = []
a0 = model.new_bool_var("a0")
e0 = model.new_bool_var("e0")
i0 = model.new_bool_var("i0")
o0 = model.new_bool_var("o0")
u0 = model.new_bool_var("u0")
model.AddExactlyOne(a0, e0, i0, o0, u0)
# ... similar "one-hot" encoding for the other free variables
n_0 = model.new_int_var(0, 2, "n_0")
model.add(n_0 == 1 + d4)
n_1 = model.new_int_var(0, 2, "n_1")
model.add(n_1 == (s3 * 2) + (y3 * 2))
n_2 = model.new_int_var(0, 3, "n_2")
model.add(n_2 == 1 + n_1)
p_0 = build_product_var(model, d4, n_2, "p_0")
n_3 = model.new_int_var(0, 3, "n_3")
model.add(n_3 == p_0)
p_1 = build_product_var(model, e2, n_0, "p_1")
p_2 = build_product_var(model, i2, n_3, "p_2")
n_4 = model.new_int_var(0, 3, "n_4")
model.add(n_4 == p_1 + p_2)
n_5 = model.new_int_var(0, 2, "n_5")
model.add(n_5 == (d4 * 2))
# ...
n_3294 = model.new_int_var(0, 45, "n_3294")
model.add(n_3294 == p_2589 + p_2590 + p_2591 + p_2592 + p_2593 + p_2594)
n_3295 = model.new_int_var(0, 250, "n_3295")
model.add(n_3295 == n_384 + n_1082 + n_1412 + n_1938 + n_2788 + n_3294)
model.add(n_3295 >= 45)
solver = cp_model.CpSolver()
status = solver.solve(model)
and finally the SMT file for Z3 (full code):
(declare-const a0 Bool)
(declare-const e0 Bool)
(declare-const i0 Bool)
(declare-const o0 Bool)
(declare-const u0 Bool)
(assert (= (+ a0 e0 i0 o0 u0) 1))
; ...
(declare-const n_0 Int)
(assert (= n_0 (+ 1 d4)))
(declare-const n_1 Int)
(assert (= n_1 (+ (* s3 2) (* y3 2))))
(declare-const n_2 Int)
(assert (= n_2 (+ 1 n_1)))
(declare-const n_3 Int)
(assert (= n_3 (+ (* d4 n_2))))
(declare-const n_4 Int)
(assert (= n_4 (+ (* e2 n_0) (* i2 n_3))))
(declare-const n_5 Int)
(assert (= n_5 (+ (* d4 2))))
; ...
(declare-const n_3294 Int)
(assert (= n_3294 (+ (* c5 n_2914) (* h5 n_2973) (* k5 n_3016) (* m5 n_3119) (* p5 n_3202) (* t5 n_3293))))
(declare-const n_3295 Int)
(assert (= n_3295 (+ n_384 n_1082 n_1412 n_1938 n_2788 n_3294)))
(assert (> n_3295 45))
(check-sat)
Here are the times I get on my M2 MacBook:
- Python brute force: 4.755s
- OR-Tools: 8.7s
- Z3: 19.5s
I would have expected brute force search to be a worst case for both solvers, and Python is hardly a fast language, so I'm wondering if there's something horribly inefficient about my encoding. A few things I've tried:
- Add
(check-sat-using (then simplify smt))for Z3: runs ~50% faster. - Use
build_product_varfrom this answer for OR-Tools. This is possibly a small speedup.
Is there something I'm missing? Or is this just a really difficult problem for these solvers?
I'm trying to establish an upper bound on a fairly elaborate equation that comes from a tree. It has a small number of of inputs that can be encoded with "one-hot" variables. All numbers are integers. All terms are sums or products of other terms and booleans. Check out this question for details.
Following a suggestion in the response to that question, I've formulated my problem using Integer Linear Programming (ILP). I've encoded the problem with Google's OR-Tools in an SMT file to feed to Microsoft's Z3. For comparison, I also implemented a brute force solution in Python.
Much to my surprise, the brute force solution in Python was by far the fastest! I would have assumed it would be by far the slowest, so what am I missing?
Here's what the Python looks like (full code):
def w(c0, c1, c2, c3, c4, c5):
a0,e0,i0,o0,u0 = c0 # exactly one of these will be True
c1,h1,k1,m1,p1,t1 = c1 # (same for all the cN variables)
a2,e2,i2,o2,u2 = c2 # ...
l3,n3,r3,s3,y3 = c3 # ...
b4,d4,f4,g4,j4,v4,w4,x4,z4 = c4 # ...
c5,h5,k5,m5,p5,t5 = c5
n_0=1+d4
n_1=(s3*2)+(y3*2)
n_2=1+n_1
n_3=(d4*n_2)
n_4=(e2*n_0)+(i2*n_3)
n_5=(d4*2)
# ...
n_3293=n_3220+n_3281+n_3292
n_3294=(c5*n_2914)+(h5*n_2973)+(k5*n_3016)+(m5*n_3119)+(p5*n_3202)+(t5*n_3293)
n_3295=n_384+n_1082+n_1412+n_1938+n_2788+n_3294
return n_3295
Here's what the OR-Tools formulation looks like (full code):
from ortools.sat.python import cp_model
model = cp_model.CpModel()
cells = []
a0 = model.new_bool_var("a0")
e0 = model.new_bool_var("e0")
i0 = model.new_bool_var("i0")
o0 = model.new_bool_var("o0")
u0 = model.new_bool_var("u0")
model.AddExactlyOne(a0, e0, i0, o0, u0)
# ... similar "one-hot" encoding for the other free variables
n_0 = model.new_int_var(0, 2, "n_0")
model.add(n_0 == 1 + d4)
n_1 = model.new_int_var(0, 2, "n_1")
model.add(n_1 == (s3 * 2) + (y3 * 2))
n_2 = model.new_int_var(0, 3, "n_2")
model.add(n_2 == 1 + n_1)
p_0 = build_product_var(model, d4, n_2, "p_0")
n_3 = model.new_int_var(0, 3, "n_3")
model.add(n_3 == p_0)
p_1 = build_product_var(model, e2, n_0, "p_1")
p_2 = build_product_var(model, i2, n_3, "p_2")
n_4 = model.new_int_var(0, 3, "n_4")
model.add(n_4 == p_1 + p_2)
n_5 = model.new_int_var(0, 2, "n_5")
model.add(n_5 == (d4 * 2))
# ...
n_3294 = model.new_int_var(0, 45, "n_3294")
model.add(n_3294 == p_2589 + p_2590 + p_2591 + p_2592 + p_2593 + p_2594)
n_3295 = model.new_int_var(0, 250, "n_3295")
model.add(n_3295 == n_384 + n_1082 + n_1412 + n_1938 + n_2788 + n_3294)
model.add(n_3295 >= 45)
solver = cp_model.CpSolver()
status = solver.solve(model)
and finally the SMT file for Z3 (full code):
(declare-const a0 Bool)
(declare-const e0 Bool)
(declare-const i0 Bool)
(declare-const o0 Bool)
(declare-const u0 Bool)
(assert (= (+ a0 e0 i0 o0 u0) 1))
; ...
(declare-const n_0 Int)
(assert (= n_0 (+ 1 d4)))
(declare-const n_1 Int)
(assert (= n_1 (+ (* s3 2) (* y3 2))))
(declare-const n_2 Int)
(assert (= n_2 (+ 1 n_1)))
(declare-const n_3 Int)
(assert (= n_3 (+ (* d4 n_2))))
(declare-const n_4 Int)
(assert (= n_4 (+ (* e2 n_0) (* i2 n_3))))
(declare-const n_5 Int)
(assert (= n_5 (+ (* d4 2))))
; ...
(declare-const n_3294 Int)
(assert (= n_3294 (+ (* c5 n_2914) (* h5 n_2973) (* k5 n_3016) (* m5 n_3119) (* p5 n_3202) (* t5 n_3293))))
(declare-const n_3295 Int)
(assert (= n_3295 (+ n_384 n_1082 n_1412 n_1938 n_2788 n_3294)))
(assert (> n_3295 45))
(check-sat)
Here are the times I get on my M2 MacBook:
- Python brute force: 4.755s
- OR-Tools: 8.7s
- Z3: 19.5s
I would have expected brute force search to be a worst case for both solvers, and Python is hardly a fast language, so I'm wondering if there's something horribly inefficient about my encoding. A few things I've tried:
- Add
(check-sat-using (then simplify smt))for Z3: runs ~50% faster. - Use
build_product_varfrom this answer for OR-Tools. This is possibly a small speedup.
Is there something I'm missing? Or is this just a really difficult problem for these solvers?
Share Improve this question asked Feb 7 at 22:14 danvkdanvk 16.9k8 gold badges78 silver badges131 bronze badges2 Answers
Reset to default 3I had a brief look into your CP-SAT code and noticed you make use of model.add_multiplication_equality quite a lot.
e.g., this is the first appearance:
n_2 = model.new_int_var(0, 3, 'n_2')
model.add(n_2 == 1+n_1)
p_0 = model.new_int_var(0, 3, 'p_0')
model.add_multiplication_equality(p_0, (d4, n_2))
While it is supported by the CP-SAT solver, it is a non-linear expression.
A formulation as a linear expression may be more favorable performance-wise.
And with …
d4 = model.new_bool_var('d4')
… it appears to be a conditional constraint with d4 as an auxiliary binary variable.
Alternatively to the big-M method to linearize such conditional constraint, CP-SAT also supports reification with only_enforce_if()
# model.add_multiplication_equality(p_0, (d4, n_2)):
model.add(p_0 == n_2).only_enforce_if(d4)
model.add(p_0 == 0).only_enforce_if(d4.Not())
# if n_2 is only used here, we can drop it as well:
# model.add(p_0 == 1+n_1).only_enforce_if(d4)
Edit:
I missed it before, just saw that you already use only_enforce_if() viabuild_product_var
Similarly, we can also linearize model.add_exactly_one()
(my experience is that it is slightly faster, whether it is significant or not, depends on the complete model)
# model.add_exactly_one(b4, d4, f4, g4, j4, v4, w4, x4, z4):
model.add(b4 + d4 + f4 + g4 + j4 + v4 + w4 + x4 + z4 == 1)
Having a ’pure’ linear formulation of your ILP problem, an LP solver might be faster than CP-SAT or Z3.
For example, with OR-Tools MathOpt you could run the ILP model against HiGHS or SCIP among others.
On the other hand, these CP-SAT constraints might be applicable for your problem as well:
- model.add_allowed_assignments()
- model.add_automaton()
Especially, the automaton constraint can speed up the solving process significantly.
Edit:
I just had a brief look at your original post of the question. These two constraints follows the same/similar idea as your point Choice Mask.And model.add_element() might also be applicable to refer to a pre-computed list of values, with a composite index based on your auxiliary variables. But I couldn’t tell whether it is beneficial for your model.
(Alternatively, these pre-computed values could also be part of the allowed_assignments)
A constructive method is expected to be faster than an exponential search.