From 88844a84aa429bb606c03fc2b4dda3ee62967fc4 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sat, 27 Sep 2025 12:17:40 +0300 Subject: [PATCH] mul-saturation wip fixup conflict explanations in mul_saturation, add parameter to enable it, add statistics --- src/math/lp/lp_params_helper.pyg | 2 +- src/math/lp/lp_settings.cpp | 1 + src/math/lp/lp_settings.h | 3 + src/math/lp/lp_types.h | 1 + src/math/lp/nla_core.cpp | 2 +- src/math/lp/nla_mul_saturate.cpp | 229 +++++++++++++++++++------------ src/math/lp/nla_mul_saturate.h | 21 +-- src/params/smt_params_helper.pyg | 1 + src/smt/theory_lra.cpp | 6 +- 9 files changed, 165 insertions(+), 101 deletions(-) diff --git a/src/math/lp/lp_params_helper.pyg b/src/math/lp/lp_params_helper.pyg index 395adf3f6..ed35e2b66 100644 --- a/src/math/lp/lp_params_helper.pyg +++ b/src/math/lp/lp_params_helper.pyg @@ -9,6 +9,6 @@ def_module_params(module_name='lp', ('dio_ignore_big_nums', BOOL, True, 'Ignore the terms with big numbers in the Diophantine handler, only relevant when dioph_eq is true'), ('dio_calls_period', UINT, 1, 'Period of calling the Diophantine handler in the final_check()'), ('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'), - ('enable_relevancy', BOOL, False, 'enabled relevancy filtering of monomials (experimental)'), + ('enable_relevancy', BOOL, False, 'enabled relevancy filtering of monomials (experimental)'), )) diff --git a/src/math/lp/lp_settings.cpp b/src/math/lp/lp_settings.cpp index 0a1a0cd54..7f3aeb2e2 100644 --- a/src/math/lp/lp_settings.cpp +++ b/src/math/lp/lp_settings.cpp @@ -34,6 +34,7 @@ void lp::lp_settings::updt_params(params_ref const& _p) { report_frequency = p.arith_rep_freq(); m_simplex_strategy = static_cast(p.arith_simplex_strategy()); m_nlsat_delay = p.arith_nl_delay(); + m_enable_stellensatz = p.arith_nl_stellensatz(); auto eps = p.arith_epsilon(); m_epsilon = rational(std::max(1, (int)(100000*eps)), 100000); m_dio = lp_p.dio(); diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h index afa73c70e..b4843e280 100644 --- a/src/math/lp/lp_settings.h +++ b/src/math/lp/lp_settings.h @@ -137,6 +137,7 @@ struct statistics { unsigned m_bounds_tightening_conflicts = 0; unsigned m_bounds_tightenings = 0; unsigned m_nla_throttled_lemmas = 0; + unsigned m_nla_stellensatz = 0; ::statistics m_st = {}; @@ -176,6 +177,7 @@ struct statistics { st.update("arith-bounds-tightening-conflicts", m_bounds_tightening_conflicts); st.update("arith-bounds-tightenings", m_bounds_tightenings); st.update("arith-nla-throttled-lemmas", m_nla_throttled_lemmas); + st.update("arith-nla-stellensatz", m_nla_stellensatz); st.copy(m_st); } }; @@ -260,6 +262,7 @@ private: bool m_dio_run_gcd = true; public: bool m_enable_relevancy = false; + bool m_enable_stellensatz = true; unsigned dio_calls_period() const { return m_dio_calls_period; } unsigned & dio_calls_period() { return m_dio_calls_period; } bool print_external_var_name() const { return m_print_external_var_name; } diff --git a/src/math/lp/lp_types.h b/src/math/lp/lp_types.h index db19e5370..9087a7a93 100644 --- a/src/math/lp/lp_types.h +++ b/src/math/lp/lp_types.h @@ -40,6 +40,7 @@ namespace lp { EQ = 0, NE = 3 }; + typedef unsigned lpvar; const lpvar null_lpvar = UINT_MAX; const constraint_index null_ci = UINT_MAX; diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index da3f60f7f..a9f8f9445 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -1333,7 +1333,7 @@ lbool core::check() { return l_false; } - if (false && no_effect()) + if (no_effect() && lp_settings().m_enable_stellensatz) ret = m_mul_saturate.saturate(); if (no_effect() && should_run_bounded_nlsat()) diff --git a/src/math/lp/nla_mul_saturate.cpp b/src/math/lp/nla_mul_saturate.cpp index e41f4491d..2dc24f1f8 100644 --- a/src/math/lp/nla_mul_saturate.cpp +++ b/src/math/lp/nla_mul_saturate.cpp @@ -1,19 +1,30 @@ /*++ Copyright (c) 2025 Microsoft Corporation - given a monic m = x * y * z ... with evaluation val(m) != val(x) * val(y) * val(z) ... + given a monic m = x * y * z ... used in a constraint that is false under the current evaluation of x,y,z saturate constraints with respect to m in other words, if a constraint contains x*y + p >= 0, - then include the constraint z >= 0 => x*y*z + z*p >= 0 + then include the constraint x*y*z + z*p >= 0 assuming current value of z is non-negative. - Check if the system with new constraints is LP feasible. - If it is not, then produce a lemma that explains the infeasibility. + Check if the system with new constraints is LP (and MIP) feasible. + If it is not, then produce a lemma that explains infeasibility. Strategy 1: The lemma is in terms of the original constraints and bounds. Strategy 2: Attempt to eliminate new monomials from the lemma by relying on Farkas multipliers. If it succeeds to eliminate new monomials we have a lemma that is a linear combination of existing variables. + The idea is that a conflict over the new system is a set of multipliers lambda, such + that lambda * A is separated from b for the constraints Axy >= b. + The coefficients in lambda are non-negative. + They correspond to variables x, y where x are variables from the input constraints + and y are new variables introduced for new monomials. + We can test if lambda allows eliminating y by taking a subset of lambda indices where + A has rows containing y_i for some fresh y_i. Replace those rows r_j by the partial sum of + of rows multiplied by lambda_j. The sum r_j1 * lambda_j1 + .. + r_jk * lambda_jk does not + contain y_i. Repeat the same process with other variables y_i'. If a sum is + generated without any y, it is a linear consequence of the new constraints but not + necessarily derivable with the old constraints. Strategy 3: The lemma uses the new constraints. --*/ @@ -38,7 +49,8 @@ namespace nla { } void mul_saturate::init_solver() { - local_solver = alloc(lp::lar_solver); + lra_solver = alloc(lp::lar_solver); + int_solver = alloc(lp::int_solver, *lra_solver); m_vars2mon.reset(); m_mon2vars.reset(); m_values.reset(); @@ -50,6 +62,7 @@ namespace nla { auto const& lra = c().lra_solver(); auto sz = lra.number_of_vars(); for (unsigned v = 0; v < sz; ++v) { + // Declare v into lra_solver unsigned w; if (m_coi.mons().contains(v)) init_monomial(v); @@ -59,24 +72,28 @@ namespace nla { auto const& t = lra.get_term(v); // Assumption: variables in coefficients are always declared before term variable. SASSERT(all_of(t, [&](auto p) { return p.j() < v; })); - w = local_solver->add_term(t.coeffs_as_vector(), v); + w = lra_solver->add_term(t.coeffs_as_vector(), v); } else - w = local_solver->add_var(v, lra.var_is_int(v)); + w = lra_solver->add_var(v, lra.var_is_int(v)); + // assert bounds on v in the new solver. VERIFY(w == v); if (lra.column_has_lower_bound(v)) { - auto lo_dep = lra.get_column_lower_bound_witness(v); auto lo_bound = lra.get_lower_bound(v); auto k = lo_bound.y > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE; - auto ci = local_solver->add_var_bound(v, k, lo_bound.x); + auto rhs = lo_bound.x; + auto dep = lra.get_column_lower_bound_witness(v); + auto ci = lra_solver->add_var_bound(v, k, rhs); + m_ci2dep.setx(ci, dep, nullptr); } if (lra.column_has_upper_bound(v)) { - auto hi_dep = lra.get_column_upper_bound_witness(v); auto hi_bound = lra.get_upper_bound(v); auto k = hi_bound.y < 0 ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE; - auto ci = local_solver->add_var_bound(v, k, hi_bound.x); - m_ci2dep.setx(ci, hi_dep, nullptr); + auto rhs = hi_bound.x; + auto dep = lra.get_column_upper_bound_witness(v); + auto ci = lra_solver->add_var_bound(v, k, rhs); + m_ci2dep.setx(ci, dep, nullptr); } } } @@ -96,81 +113,106 @@ namespace nla { void mul_saturate::add_lemma(lp::explanation const& ex1) { auto& lra = c().lra_solver(); lp::explanation ex2; - for (auto p : ex1) { - lp::constraint_index src_ci; - if (!m_new_mul_constraints.find(p.ci(), src_ci)) - src_ci = p.ci(); - auto dep = m_ci2dep.get(src_ci, nullptr); - local_solver->push_explanation(dep, ex2); - } - for (auto [v, sign, dep] : m_var_signs) { - if (!dep) { - verbose_stream() << "TODO explain non-implied bound\n"; - continue; - } - local_solver->push_explanation(dep, ex2); - } lemma_builder new_lemma(c(), "stellensatz"); + for (auto p : ex1) { + auto dep = m_ci2dep.get(p.ci(), nullptr); + lra_solver->push_explanation(dep, ex2); + if (!m_new_mul_constraints.contains(p.ci())) + continue; + + auto const& bounds = m_new_mul_constraints[p.ci()]; + for (auto const& b : bounds) { + if (std::holds_alternative(b)) { + auto dep = *std::get_if(&b); + lra_solver->push_explanation(dep, ex2); + } + else { + auto const &[v, k, rhs] = *std::get_if(&b); + new_lemma |= ineq(v, negate(k), rhs); + } + } + } + new_lemma &= ex2; - IF_VERBOSE(1, verbose_stream() << "unsat \n" << new_lemma << "\n"); + IF_VERBOSE(5, verbose_stream() << "unsat \n" << new_lemma << "\n"); TRACE(arith, tout << "unsat\n" << new_lemma << "\n"); + c().lra_solver().settings().stats().m_nla_stellensatz++; } + lbool mul_saturate::solve(lp::explanation& ex) { - for (auto const& [new_ci, old_ci] : m_new_mul_constraints) - local_solver->activate(new_ci); - auto st = local_solver->solve(); - lbool r = l_undef; + lbool r = solve_lra(ex); + if (r != l_true) + return r; + r = solve_lia(ex); + if (r != l_true) + return r; + // if (r == l_true) check if solution satisfies constraints + // variables outside the slice have values from the outer solver. + return l_undef; + } + + lbool mul_saturate::solve_lra(lp::explanation& ex) { + auto st = lra_solver->solve(); if (st == lp::lp_status::INFEASIBLE) { - local_solver->get_infeasibility_explanation(ex); - r = l_false; + lra_solver->get_infeasibility_explanation(ex); + return l_false; + } + else if (lra_solver->is_feasible()) + return l_true; + else + return l_undef; + } + + lbool mul_saturate::solve_lia(lp::explanation& ex) { + switch (int_solver->check(&ex)) { + case lp::lia_move::sat: + return l_true; + case lp::lia_move::conflict: + return l_false; + default: // TODO: an option is to perform (bounded) search here to get an LIA verdict. + return l_undef; } - if (st == lp::lp_status::OPTIMAL || st == lp::lp_status::FEASIBLE) { - // TODO: check model just in case it got lucky. - IF_VERBOSE(1, verbose_stream() << "saturation is LP feasible, maybe it is a model of the NLA problem\n"); - } - IF_VERBOSE(0, display(verbose_stream())); - return r; + return l_undef; } // record new monomials that are created and recursively down-saturate with respect to these. // this is a simplistic pass void mul_saturate::add_multiply_constraints() { m_new_mul_constraints.reset(); - m_seen_vars.reset(); - m_var_signs.reset(); m_to_refine.reset(); - vector> var2cs; + vector> var2cs; - for (auto ci : local_solver->constraints().indices()) { - auto const& con = local_solver->constraints()[ci]; + // current approach: only resolve against var2cs, which is initialized + // with monomials in the input. + + for (auto ci : lra_solver->constraints().indices()) { + auto const& con = lra_solver->constraints()[ci]; for (auto [coeff, v] : con.coeffs()) { if (v >= var2cs.size()) var2cs.resize(v + 1); var2cs[v].push_back(ci); } - // insert monomials to be refined insert_monomials_from_constraint(ci); } for (unsigned it = 0; it < m_to_refine.size(); ++it) { auto j = m_to_refine[it]; - verbose_stream() << "refining " << j << " := " << m_mon2vars[j] << "\n"; auto vars = m_mon2vars[j]; for (auto v : vars) { if (v >= var2cs.size()) continue; auto cs = var2cs[v]; for (auto ci : cs) { - for (auto [coeff, u] : local_solver->constraints()[ci].coeffs()) { + for (auto [coeff, u] : lra_solver->constraints()[ci].coeffs()) { if (u == v) add_multiply_constraint(ci, j, v); } } } } - IF_VERBOSE(0, + IF_VERBOSE(5, c().lra_solver().display(verbose_stream() << "original\n"); c().display(verbose_stream()); display(verbose_stream() << "saturated\n")); @@ -178,52 +220,60 @@ namespace nla { // multiply by remaining vars void mul_saturate::add_multiply_constraint(lp::constraint_index old_ci, lp::lpvar mi, lpvar x) { - lp::lar_base_constraint const& con = local_solver->constraints()[old_ci]; + lp::lar_base_constraint const& con = lra_solver->constraints()[old_ci]; auto &lra = c().lra_solver(); auto const& lhs = con.coeffs(); auto const& rhs = con.rhs(); auto k = con.kind(); if (k == lp::lconstraint_kind::NE || k == lp::lconstraint_kind::EQ) return; // not supported - auto sign = false; + auto sign = false, strict = true; svector vars; bool first = true; - for (auto v : c().emon(mi).vars()) { + for (auto v : m_mon2vars[mi]) { if (v != x || !first) vars.push_back(v); else first = false; } - // compute sign of vars + SASSERT(!first); // v was a member and was removed + + vector bounds; + // compute bounds constraints and sign of vars + + auto add_bound = [&](lpvar v, lp::lconstraint_kind k, int r) { + bound b(v, k, rational(r)); + bounds.push_back(b); + }; for (auto v : vars) { - if (m_seen_vars.contains(v)) - continue; - m_seen_vars.insert(v); // retrieve bounds of v // if v has positive lower bound add as positive // if v has negative upper bound add as negative - // otherwise, soft-fail (for now unsound) - // proper signs of variables from old tableau should be extracted using lra_solver() - // instead of local_solver. - // TODO is to also add case where lower or upper bound is zero and then the sign - // of the inequality is relaxed if it is strict. + // otherwise look at the current value of v and add bounds assumption based on current sign. if (lra.number_of_vars() > v && lra.column_has_lower_bound(v) && lra.get_lower_bound(v).is_pos()) { - m_var_signs.push_back({v, false, lra.get_column_lower_bound_witness(v)}); + bounds.push_back(lra.get_column_lower_bound_witness(v)); } else if (lra.number_of_vars() > v && lra.column_has_upper_bound(v) && lra.get_upper_bound(v).is_neg()) { - m_var_signs.push_back({v, true, lra.get_column_upper_bound_witness(v)}); + bounds.push_back(lra.get_column_upper_bound_witness(v)); sign = !sign; } else if (m_values[v].is_neg()) { - m_var_signs.push_back({v, true, nullptr}); + if (lra.var_is_int(v)) + add_bound(v, lp::lconstraint_kind::LE, -1); + else + add_bound(v, lp::lconstraint_kind::LT, 0); sign = !sign; } else if (m_values[v].is_pos()) { - m_var_signs.push_back({v, false, nullptr}); + if (lra.var_is_int(v)) + add_bound(v, lp::lconstraint_kind::GE, 1); + else + add_bound(v, lp::lconstraint_kind::GT, 0); } else { - IF_VERBOSE(0, verbose_stream() << "not separated from 0\n"); - return; + SASSERT(m_values[v] == 0); + strict = false; + add_bound(v, lp::lconstraint_kind::GE, 0); } } lp::lar_term new_lhs; @@ -246,34 +296,35 @@ namespace nla { new_rhs = 0; } - if (sign) { + if (sign) + k = swap_side(k); + + if (!strict) { switch (k) { - case lp::lconstraint_kind::LE: + case lp::lconstraint_kind::GT: k = lp::lconstraint_kind::GE; break; case lp::lconstraint_kind::LT: - k = lp::lconstraint_kind::GT; - break; - case lp::lconstraint_kind::GE: k = lp::lconstraint_kind::LE; break; - case lp::lconstraint_kind::GT: - k = lp::lconstraint_kind::LT; - break; - default: - break; } } - display_constraint(verbose_stream(), old_ci) << " -> "; - display_constraint(verbose_stream(), new_lhs.coeffs_as_vector(), k, new_rhs) << "\n"; - auto ti = local_solver->add_term(new_lhs.coeffs_as_vector(), local_solver->number_of_vars()); - auto new_ci = local_solver->add_var_bound(ti, k, new_rhs); + + auto ti = lra_solver->add_term(new_lhs.coeffs_as_vector(), lra_solver->number_of_vars()); + auto new_ci = lra_solver->add_var_bound(ti, k, new_rhs); + IF_VERBOSE(4, display_constraint(verbose_stream(), old_ci) << " -> "; + display_constraint(verbose_stream(), new_lhs.coeffs_as_vector(), k, new_rhs) << "\n"); insert_monomials_from_constraint(new_ci); m_values.push_back(term_value); SASSERT(m_values.size() - 1 == ti); - m_new_mul_constraints.insert(new_ci, old_ci); + m_new_mul_constraints.insert(new_ci, bounds); + m_ci2dep.setx(new_ci, m_ci2dep.get(old_ci, nullptr), nullptr); } + // Insert a (new) monomial representing product of vars. + // if the length of vars is 1 then the new monomial is vars[0] + // if the same monomial was previous defined, return the previously defined monomial. + // otherwise create a fresh variable representing the monomial. lpvar mul_saturate::add_monomial(svector const& vars) { lpvar v; if (vars.size() == 1) @@ -300,23 +351,25 @@ namespace nla { } lpvar mul_saturate::add_var(bool is_int) { - auto v = local_solver->number_of_vars(); - auto w = local_solver->add_var(v, is_int); - VERIFY(v == w); + auto v = lra_solver->number_of_vars(); + auto w = lra_solver->add_var(v, is_int); + SASSERT(v == w); return w; } + // if a constraint is false, then insert _all_ monomials from that constraint. + // other approach: use a lex ordering on monomials and insert the lex leading monomial. void mul_saturate::insert_monomials_from_constraint(lp::constraint_index ci) { if (constraint_is_true(ci)) return; - auto const& con = local_solver->constraints()[ci]; + auto const& con = lra_solver->constraints()[ci]; for (auto [coeff, v] : con.coeffs()) if (c().is_monic_var(v)) m_to_refine.insert(v); } bool mul_saturate::constraint_is_true(lp::constraint_index ci) { - auto const& con = local_solver->constraints()[ci]; + auto const& con = lra_solver->constraints()[ci]; auto lhs = -con.rhs(); for (auto const& [coeff, v] : con.coeffs()) lhs += coeff * m_values[v]; @@ -340,7 +393,7 @@ namespace nla { } std::ostream& mul_saturate::display(std::ostream& out) const { - local_solver->display(out); + lra_solver->display(out); for (auto const& [vars, v] : m_vars2mon) { out << "j" << v << " := "; display_product(out, vars); @@ -362,7 +415,7 @@ namespace nla { } std::ostream& mul_saturate::display_constraint(std::ostream& out, lp::constraint_index ci) const { - auto const& con = local_solver->constraints()[ci]; + auto const& con = lra_solver->constraints()[ci]; return display_constraint(out, con.coeffs(), con.kind(), con.rhs()); } diff --git a/src/math/lp/nla_mul_saturate.h b/src/math/lp/nla_mul_saturate.h index 0c5e09e1f..247fc8381 100644 --- a/src/math/lp/nla_mul_saturate.h +++ b/src/math/lp/nla_mul_saturate.h @@ -5,6 +5,8 @@ #pragma once #include "math/lp/nla_coi.h" +#include "math/lp/int_solver.h" +#include namespace nla { @@ -12,18 +14,19 @@ namespace nla { class lar_solver; class mul_saturate : common { - struct var_sign { + + struct bound { lpvar v = lp::null_lpvar; - bool is_neg = false; - u_dependency* dep = nullptr; + lp::lconstraint_kind k; + rational rhs; }; + using bound_justification = std::variant; + coi m_coi; - // source of multiplication constraint - u_map m_new_mul_constraints; - svector m_var_signs; - tracked_uint_set m_seen_vars; + u_map> m_new_mul_constraints; indexed_uint_set m_to_refine; - scoped_ptr local_solver; + scoped_ptr lra_solver; + scoped_ptr int_solver; ptr_vector m_ci2dep; vector m_values; struct eq { @@ -51,6 +54,8 @@ namespace nla { // solving lbool solve(lp::explanation& ex); + lbool solve_lra(lp::explanation &ex); + lbool solve_lia(lp::explanation &ex); // lemmas void add_lemma(lp::explanation const& ex); diff --git a/src/params/smt_params_helper.pyg b/src/params/smt_params_helper.pyg index 487772c81..58d08ce8a 100644 --- a/src/params/smt_params_helper.pyg +++ b/src/params/smt_params_helper.pyg @@ -88,6 +88,7 @@ def_module_params(module_name='smt', ('arith.nl.propagate_linear_monomials', BOOL, True, 'propagate linear monomials'), ('arith.nl.optimize_bounds', BOOL, True, 'enable bounds optimization'), ('arith.nl.cross_nested', BOOL, True, 'enable cross-nested consistency checking'), + ('arith.nl.stellensatz', BOOL, False, 'enable stellensatz saturation'), ('arith.nl.log', BOOL, False, 'Log lemmas sent to nra solver'), ('arith.propagate_eqs', BOOL, True, 'propagate (cheap) equalities'), ('arith.epsilon', DOUBLE, 1.0, 'initial value of epsilon used for model generation of infinitesimals'), diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 6876d363c..1d94fbcd2 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -3467,15 +3467,15 @@ public: void set_conflict_or_lemma(literal_vector const& core, bool is_conflict) { reset_evidence(); - for (literal lit : core) { + for (literal lit : core) m_core.push_back(lit); - } + // lp().shrink_explanation_to_minimum(m_explanation); // todo, enable when perf is fixed ++m_num_conflicts; ++m_stats.m_conflicts; for (auto ev : m_explanation) set_evidence(ev.ci(), m_core, m_eqs); - if (m_eqs.empty() && all_of(m_core, [&](literal l) { return ctx().get_assignment(l) == l_false; })) + if (all_of(m_core, [&](literal l) { return ctx().get_assignment(l) == l_false; })) is_conflict = true; TRACE(arith_conflict, tout << "@" << ctx().get_scope_level() << (is_conflict ? " conflict":" lemma");