From bed532442b3378de9ff1455eb322b690eef3fb95 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Sat, 16 May 2026 07:31:03 -0700 Subject: [PATCH] lll-cube: restore C(B) monotonicity, hoist invariants, cap candidate q Three combined fixes to the lexicographic basis-reduction heuristic. 1. Monotonicity guard. The lexicographic scoring in reduce_pair can accept moves with cost > cost(B=I) if they reduce the bound-collapse excess. Across passes this lets the basis drift to C(B) > C(I), which then bloats compute_deltas and triggers bail_tighten or wastes work in find_feasible_solution. We maintain an unweighted running C(B) and compare against the C(I) baseline taken right after build_matrix; if the basis search did not strictly improve, operator() bails as lia_move::undef. Restores the pre-lex 'never worse than plain int_cube' invariant noted in the header. 2. Candidate magnitude cap. The new candidate set includes values of the form floor((a +/- (width - base)) / b), where 'width' is the raw bound interval. Wide-box columns (e.g. width = 1e6) produced q values far outside the breakpoint range, seeding basis growth that overflowed too_big() downstream. Add an explicit 64-bit bitsize ceiling on candidates, well below LLL_CUBE_BITSIZE_LIMIT. 3. Hoist invariants out of reduce_pair. Eligibility, width, and integral_unit are computed once per operator() call (init_pair_- invariants). Row 1-norms are cached in m_row_sum_H/m_row_sum_B and maintained incrementally on each accepted move; the inner base = sum_{c != j} |M[r][c]| computation collapses from O(n) per risk row per pair to O(1) (row_sum - |entry|). Per-pair cost in the risk-row build drops from O((m+n)*n) to O(m+n); at the limits (m=75, n=150, 8 passes, ~27 candidates) this is ~6-7x fewer mpq ops per compute_basis call, with the dominant savings being the per-row 'mpq base(0); base += abs(...)' GMP allocations. Soundness is unchanged: the moves applied are still unimodular column ops, and the rounding path still uses B^{-1} and B as before. Signed-off-by: Lev Nachmanson Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- src/math/lp/int_cube_lll.cpp | 140 +++++++++++++++++++++++++++-------- src/math/lp/int_cube_lll.h | 22 ++++++ 2 files changed, 130 insertions(+), 32 deletions(-) diff --git a/src/math/lp/int_cube_lll.cpp b/src/math/lp/int_cube_lll.cpp index c49d0a677..5328befd5 100644 --- a/src/math/lp/int_cube_lll.cpp +++ b/src/math/lp/int_cube_lll.cpp @@ -26,6 +26,7 @@ namespace lp { static const unsigned LLL_CUBE_LIMIT_ROWS = 75; static const unsigned LLL_CUBE_LIMIT_COLS = 150; static const unsigned LLL_CUBE_BITSIZE_LIMIT = 4096; + static const unsigned LLL_CUBE_CANDIDATE_BITSIZE = 64; static const unsigned LLL_CUBE_MAX_PASSES = 8; int_cube_lll::int_cube_lll(int_solver& lia) : lia(lia), lra(lia.lra) {} @@ -118,9 +119,70 @@ namespace lp { m_B_inv[i][i] = mpq(1); } compute_col_weights(); + init_pair_invariants(); return true; } + // Cache per-row metadata, row 1-norms, and the identity-baseline cube + // cost C(I). Invariants of lra/bounds for the duration of operator(), + // so we pay this once instead of per (j, k) reduce_pair call. + void int_cube_lll::init_pair_invariants() { + unsigned m = m_H.size(); + unsigned n = m_B.size(); + m_term_meta.reset(); + m_term_meta.resize(m); + for (unsigned r = 0; r < m; ++r) { + unsigned var = m_term_js[r]; + row_meta& meta = m_term_meta[r]; + meta.eligible = lra.column_has_lower_bound(var) && lra.column_has_upper_bound(var); + if (meta.eligible) { + meta.width = lra.get_upper_bound(var).x - lra.get_lower_bound(var).x; + meta.integral_unit = column_bounds_are_integral(var); + } + else { + meta.width = mpq(0); + meta.integral_unit = false; + } + } + m_col_meta.reset(); + m_col_meta.resize(n); + for (unsigned i = 0; i < n; ++i) { + unsigned var = m_J[i]; + row_meta& meta = m_col_meta[i]; + meta.eligible = lra.column_has_lower_bound(var) && lra.column_has_upper_bound(var); + if (meta.eligible) { + meta.width = lra.get_upper_bound(var).x - lra.get_lower_bound(var).x; + meta.integral_unit = column_bounds_are_integral(var); + } + else { + meta.width = mpq(0); + meta.integral_unit = false; + } + } + m_row_sum_H.reset(); + m_row_sum_H.resize(m, mpq(0)); + for (unsigned r = 0; r < m; ++r) { + mpq s(0); + for (unsigned c = 0; c < n; ++c) + s += abs(m_H[r][c]); + m_row_sum_H[r] = s; + } + m_row_sum_B.reset(); + m_row_sum_B.resize(n, mpq(0)); + for (unsigned i = 0; i < n; ++i) { + mpq s(0); + for (unsigned c = 0; c < n; ++c) + s += abs(m_B[i][c]); + m_row_sum_B[i] = s; + } + m_baseline_cost = mpq(0); + for (unsigned r = 0; r < m; ++r) + m_baseline_cost += m_row_sum_H[r]; + for (unsigned i = 0; i < n; ++i) + m_baseline_cost += m_row_sum_B[i]; + m_total_cost = m_baseline_cost; + } + // Weight per column for the cube-cost function used by the greedy // basis reduction. A tight-box column (small ub - lb) is expensive // to grow because tightening by delta_col(j) >= (ub - lb)/2 would @@ -221,45 +283,37 @@ namespace lp { mpq qf = floor(median); mpq qc = qf + mpq(1); + // Risk rows are built in O(m + n) using cached metadata and the + // maintained row 1-norms. base(r, j) = ||row_r||_1 - |row_r[j]| is + // the partial row sum that would survive the column-j update. struct risk_row { - unsigned var; + const row_meta* meta; mpq base; mpq a; mpq b; - mpq width; - bool integral_unit; }; vector risk_rows; + risk_rows.reserve(m + n); for (unsigned r = 0; r < m; ++r) { - unsigned var = m_term_js[r]; - if (!lra.column_has_lower_bound(var) || !lra.column_has_upper_bound(var)) + const row_meta& meta = m_term_meta[r]; + if (!meta.eligible) continue; - mpq base(0); - for (unsigned c = 0; c < n; ++c) - if (c != j) - base += abs(m_H[r][c]); - risk_rows.push_back({var, base, m_H[r][j], m_H[r][k], - lra.get_upper_bound(var).x - lra.get_lower_bound(var).x, - column_bounds_are_integral(var)}); + mpq base = m_row_sum_H[r] - abs(m_H[r][j]); + risk_rows.push_back({&meta, base, m_H[r][j], m_H[r][k]}); } for (unsigned i = 0; i < n; ++i) { - unsigned var = m_J[i]; - if (!lra.column_has_lower_bound(var) || !lra.column_has_upper_bound(var)) + const row_meta& meta = m_col_meta[i]; + if (!meta.eligible) continue; - mpq base(0); - for (unsigned c = 0; c < n; ++c) - if (c != j) - base += abs(m_B[i][c]); - risk_rows.push_back({var, base, m_B[i][j], m_B[i][k], - lra.get_upper_bound(var).x - lra.get_lower_bound(var).x, - column_bounds_are_integral(var)}); + mpq base = m_row_sum_B[i] - abs(m_B[i][j]); + risk_rows.push_back({&meta, base, m_B[i][j], m_B[i][k]}); } auto row_excess = [](const risk_row& r, const mpq& q) { mpq norm = r.base + abs(r.a - q * r.b); - if (norm.is_zero() || (norm == mpq(1) && r.integral_unit)) + if (norm.is_zero() || (norm == mpq(1) && r.meta->integral_unit)) return mpq(0); - return norm > r.width ? norm - r.width : mpq(0); + return norm > r.meta->width ? norm - r.meta->width : mpq(0); }; struct score { @@ -288,8 +342,11 @@ namespace lp { }; vector candidates; + // Cap candidate magnitude well below LLL_CUBE_BITSIZE_LIMIT so that + // floor((a +/- width)/b) values for wide-box columns don't seed + // basis growth and overflow downstream. auto add_candidate = [&](const mpq& q) { - if (too_big(q)) + if (q.bitsize() > LLL_CUBE_CANDIDATE_BITSIZE) return; for (auto const& c : candidates) if (c == q) @@ -322,7 +379,7 @@ namespace lp { mpq f = floor(center); add_candidate(f); add_candidate(f + mpq(1)); - mpq allowance = r.width - r.base; + mpq allowance = r.meta->width - r.base; if (allowance >= mpq(0)) { mpq lo = (r.a - allowance) / r.b; mpq hi = (r.a + allowance) / r.b; @@ -350,19 +407,29 @@ namespace lp { } // Apply column op col_j <- col_j - q * col_k on H and B, and the - // matching row op row_k <- row_k + q * row_j on B_inv. + // matching row op row_k <- row_k + q * row_j on B_inv. Update the + // cached row 1-norms and running total cost in lock-step. const mpq& q = best_q; + mpq delta_total(0); for (unsigned r = 0; r < m; ++r) { - mpq v = m_H[r][j] - q * m_H[r][k]; - if (too_big(v)) + mpq old_a = m_H[r][j]; + mpq new_a = old_a - q * m_H[r][k]; + if (too_big(new_a)) return false; - m_H[r][j] = v; + m_H[r][j] = new_a; + mpq d = abs(new_a) - abs(old_a); + m_row_sum_H[r] += d; + delta_total += d; } for (unsigned i = 0; i < n; ++i) { - mpq v = m_B[i][j] - q * m_B[i][k]; - if (too_big(v)) + mpq old_a = m_B[i][j]; + mpq new_a = old_a - q * m_B[i][k]; + if (too_big(new_a)) return false; - m_B[i][j] = v; + m_B[i][j] = new_a; + mpq d = abs(new_a) - abs(old_a); + m_row_sum_B[i] += d; + delta_total += d; } for (unsigned c = 0; c < n; ++c) { mpq v = m_B_inv[k][c] + q * m_B_inv[j][c]; @@ -370,6 +437,7 @@ namespace lp { return false; m_B_inv[k][c] = v; } + m_total_cost += delta_total; improved = true; return true; } @@ -549,6 +617,14 @@ namespace lp { lia.settings().stats().m_lll_cube_bail_basis++; return lia_move::undef; } + // Monotonicity guard: bail if the computed basis did not strictly + // improve the unweighted cube cost C(B) below the identity baseline + // C(I) = ||A||_1 + n. This restores the "never worse than plain + // int_cube" guarantee that pre-lexicographic-scoring code provided. + if (m_total_cost >= m_baseline_cost) { + lia.settings().stats().m_lll_cube_bail_basis++; + return lia_move::undef; + } compute_deltas(); lra.push(); diff --git a/src/math/lp/int_cube_lll.h b/src/math/lp/int_cube_lll.h index e102a2c27..8f1ebb71b 100644 --- a/src/math/lp/int_cube_lll.h +++ b/src/math/lp/int_cube_lll.h @@ -28,6 +28,10 @@ Abstract: of bounded rows/columns whose later delta-tightening would collapse their box, then minimize C(B). This keeps the greedy from accepting a cost-improving basis move that merely sets up a bail-tighten failure. + Soundness: after the basis search completes we compare the unweighted + C(B) against the identity baseline C(I) = ||A||_1 + n; if the basis + failed to strictly improve it we bail back to plain int_cube semantics. + The heuristic is therefore never worse than the plain int_cube. This addresses the regression of int_cube_hnf, whose triangulation can blow up the column-delta term ||B||_1: in 153 random instances @@ -62,6 +66,17 @@ namespace lp { class int_solver& lia; class lar_solver& lra; + // Per-row/column metadata cached for the duration of one operator() + // call. All fields are invariants of the LP bounds and lattice + // structure (lra is read-only during compute_basis), so they are + // computed once in init_pair_invariants() instead of being rebuilt + // inside every reduce_pair(j, k) call. + struct row_meta { + bool eligible; // column_has_lower_bound && column_has_upper_bound + mpq width; // ub - lb (when eligible) + bool integral_unit; // column_bounds_are_integral(var) + }; + vector m_J; vector m_col_to_J; vector m_term_js; @@ -69,6 +84,12 @@ namespace lp { vector> m_B; // n x n unimodular vector> m_B_inv; // B^{-1} = product of inverse elementaries vector m_col_w; // per-column weight in the cube cost + vector m_term_meta; // per term row of H + vector m_col_meta; // per J column (row of B) + vector m_row_sum_H; // ||row_r(H)||_1, maintained incrementally + vector m_row_sum_B; // ||row_i(B)||_1, maintained incrementally + mpq m_baseline_cost; // unweighted C(B) at B = I + mpq m_total_cost; // current unweighted C(B) vector m_term_delta; vector m_col_delta; vector m_saved_x_J; @@ -81,6 +102,7 @@ namespace lp { bool collect_J_and_terms(); bool build_matrix(); void compute_col_weights(); + void init_pair_invariants(); bool compute_basis(); bool reduce_pair(unsigned j, unsigned k, bool& improved); bool too_big(const mpq& v) const;