diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ index 11ad8ca68..b00d592a3 100644 --- a/docs/paper/reductions.typ +++ b/docs/paper/reductions.typ @@ -4336,6 +4336,65 @@ Each reduction is presented as a *Rule* (with linked problem names and overhead _Solution extraction._ Convert binary to spins: $s_i = 2x_i - 1$, i.e.\ $x_i = 1 arrow.r s_i = +1$, $x_i = 0 arrow.r s_i = -1$. ] +#let cvp_qubo = load-example("ClosestVectorProblem", "QUBO") +#let cvp_qubo_sol = cvp_qubo.solutions.at(0) +#{ + let basis = cvp_qubo.source.instance.basis + let bounds = cvp_qubo.source.instance.bounds + let target = cvp_qubo.source.instance.target + let offsets = cvp_qubo_sol.source_config + let coords = offsets.enumerate().map(((i, off)) => off + bounds.at(i).lower) + let matrix = cvp_qubo.target.instance.matrix + let bits = cvp_qubo_sol.target_config + let lo = bounds.map(b => b.lower) + let anchor = range(target.len()).map(d => lo.enumerate().fold(0.0, (acc, (i, x)) => acc + x * basis.at(i).at(d))) + let constant = range(target.len()).fold(0.0, (acc, d) => acc + calc.pow(anchor.at(d) - target.at(d), 2)) + let qubo-value = range(bits.len()).fold(0.0, (acc, i) => acc + if bits.at(i) == 0 { 0.0 } else { + range(bits.len() - i).fold(0.0, (row-acc, delta) => row-acc + if bits.at(i + delta) == 0 { 0.0 } else { matrix.at(i).at(i + delta) }) + }) + let fmt-vec(v) = $paren.l #v.map(e => str(e)).join(", ") paren.r^top$ + let rounded-constant = calc.round(constant, digits: 2) + let rounded-qubo = calc.round(qubo-value, digits: 1) + let rounded-distance-sq = calc.round(qubo-value + constant, digits: 2) + [ + #reduction-rule("ClosestVectorProblem", "QUBO", + example: true, + example-caption: [2D bounded CVP with two 3-bit exact-range encodings], + extra: [ + *Step 1 -- Source instance.* The canonical CVP example uses basis columns $bold(b)_1 = #fmt-vec(basis.at(0))$ and $bold(b)_2 = #fmt-vec(basis.at(1))$, target $bold(t) = #fmt-vec(target)$, and bounds $x_1, x_2 in [#bounds.at(0).lower, #bounds.at(0).upper]$. + + *Step 2 -- Exact bounded encoding.* Each variable has #bounds.at(0).upper - bounds.at(0).lower + 1 admissible values, so the implementation uses the capped binary basis $(1, 2, 3)$ rather than $(1, 2, 4)$: the first two bits are powers of two, and the last weight is capped so every bit pattern reconstructs an offset in ${0, dots, 6}$. Thus + $ x_1 = #bounds.at(0).lower + z_0 + 2 z_1 + 3 z_2, quad x_2 = #bounds.at(1).lower + z_3 + 2 z_4 + 3 z_5 $ + giving #cvp_qubo.target.instance.num_vars QUBO variables in total. + + *Step 3 -- Build the QUBO.* For this instance, $G = A^top A = ((4, 2), (2, 5))$ and $h = A^top bold(t) = (5.6, 5.8)^top$. Expanding the shifted quadratic form yields the exported upper-triangular matrix with representative entries $Q_(0,0) = #matrix.at(0).at(0)$, $Q_(0,1) = #matrix.at(0).at(1)$, $Q_(0,2) = #matrix.at(0).at(2)$, $Q_(2,5) = #matrix.at(2).at(5)$, and $Q_(5,5) = #matrix.at(5).at(5)$. + + *Step 4 -- Verify a solution.* The fixture stores the canonical witness $bold(z) = (#bits.map(str).join(", "))$, which extracts to source offsets $bold(c) = (#offsets.map(str).join(", "))$ and actual lattice coordinates $bold(x) = (#coords.map(str).join(", "))$. The QUBO value is $bold(z)^top Q bold(z) = #rounded-qubo$; adding back the dropped constant #rounded-constant yields the original squared distance #(rounded-distance-sq), so the extracted point is the closest lattice vector #sym.checkmark. + + *Multiplicity.* Offset $3$ has two bit encodings ($(0, 0, 1)$ and $(1, 1, 0)$), so the fixture stores one canonical witness even though the QUBO has multiple optimal binary assignments representing the same CVP solution. + ], + )[ + A bounded Closest Vector Problem instance already supplies a finite integer box $x_i in [ell_i, u_i]$ for each coefficient. Following the direct quadratic-form reduction of Canale, Qureshi, and Viola @canale2023qubo, encoding each offset $c_i = x_i - ell_i$ with an exact in-range binary basis turns the squared-distance objective into an unconstrained quadratic over binary variables. Unlike penalty-method encodings, no auxiliary feasibility penalty is needed: every bit pattern decodes to a legal coefficient vector by construction. + ][ + _Construction._ Let $A in ZZ^(m times n)$ be the basis matrix with columns $bold(a)_1, dots, bold(a)_n$, let $bold(t) in RR^m$ be the target, and let $x_i in [ell_i, u_i]$ with range $r_i = u_i - ell_i$. Define $L_i = ceil(log_2(r_i + 1))$ when $r_i > 0$ and omit bits when $r_i = 0$. For each variable, introduce binary variables $z_(i,0), dots, z_(i,L_i-1)$ with exact-range weights + $ w_(i,p) = 2^p quad (0 <= p < L_i - 1), quad w_(i,L_i-1) = r_i + 1 - 2^(L_i - 1) $ + so that every bit vector represents an offset in ${0, dots, r_i}$. Then + $ x_i = ell_i + sum_(p=0)^(L_i-1) w_(i,p) z_(i,p) $ + and the total number of QUBO variables is $N = sum_i L_i$, exactly the exported overhead `num_vars = num_encoding_bits`. + + Let $G = A^top A$ and $h = A^top bold(t)$. Writing $bold(x) = bold(ell) + B bold(z)$ for the encoding matrix $B in RR^(n times N)$ gives + $ norm(A bold(x) - bold(t))_2^2 = bold(z)^top (B^top G B) bold(z) + 2 bold(z)^top B^top (G bold(ell) - h) + "const" $ + where the constant $norm(A bold(ell) - bold(t))_2^2$ is dropped. Therefore the QUBO coefficients are + $ Q_(u,u) = (B^top G B)_(u,u) + 2 (B^top (G bold(ell) - h))_u, quad Q_(u,v) = 2 (B^top G B)_(u,v) quad (u < v) $ + using the usual upper-triangular convention. + + _Correctness._ ($arrow.r.double$) Every binary vector $bold(z) in {0,1}^N$ decodes to a coefficient vector $bold(x)$ inside the prescribed bounds because each exact-range basis reaches only offsets in ${0, dots, r_i}$. Substituting this decoding into the CVP objective yields $bold(z)^top Q bold(z) + "const"$, so any QUBO minimizer maps to a bounded CVP minimizer. ($arrow.l.double$) Every bounded CVP solution $bold(x)$ has at least one bit encoding for each coordinate offset, hence at least one binary vector $bold(z)$ with the same objective value up to the dropped constant. Thus the minimizers correspond exactly, although several binary witnesses may decode to the same CVP solution. + + _Solution extraction._ For each source variable, sum its selected encoding weights to recover the source configuration offset $c_i = x_i - ell_i$. This is exactly the configuration format expected by the `ClosestVectorProblem` model. + ] + ] +} + == Penalty-Method QUBO Reductions The _penalty method_ @glover2019 @lucas2014 converts a constrained optimization problem into an unconstrained QUBO by adding quadratic penalty terms. Given an objective $"obj"(bold(x))$ to minimize and constraints $g_k (bold(x)) = 0$, construct: @@ -5250,6 +5309,7 @@ The following table shows concrete variable overhead for example instances, take (source: "SpinGlass", target: "MaxCut"), (source: "SpinGlass", target: "QUBO"), (source: "QUBO", target: "SpinGlass"), + (source: "ClosestVectorProblem", target: "QUBO"), (source: "KColoring", target: "QUBO"), (source: "MaximumSetPacking", target: "QUBO"), ( diff --git a/docs/paper/references.bib b/docs/paper/references.bib index 54f96f390..9b87db469 100644 --- a/docs/paper/references.bib +++ b/docs/paper/references.bib @@ -507,6 +507,15 @@ @article{pan2025 archivePrefix = {arXiv} } +@article{canale2023qubo, + author = {Eduardo Canale and Claudio Qureshi and Alfredo Viola}, + title = {Qubo model for the Closest Vector Problem}, + journal = {arXiv preprint}, + year = {2023}, + eprint = {2304.03616}, + archivePrefix = {arXiv} +} + @article{goemans1995, author = {Michel X. Goemans and David P. Williamson}, title = {Improved Approximation Algorithms for Maximum Cut and Satisfiability Problems Using Semidefinite Programming}, @@ -793,17 +802,6 @@ @article{cygan2014 doi = {10.1137/140990255} } -@article{dreyfuswagner1971, - author = {Stuart E. Dreyfus and Robert A. Wagner}, - title = {The Steiner Problem in Graphs}, - journal = {Networks}, - volume = {1}, - number = {3}, - pages = {195--207}, - year = {1971}, - doi = {10.1002/net.3230010302} -} - @inproceedings{bjorklund2007, author = {Andreas Bj\"{o}rklund and Thore Husfeldt and Petteri Kaski and Mikko Koivisto}, title = {Fourier Meets M\"{o}bius: Fast Subset Convolution}, diff --git a/src/models/algebraic/closest_vector_problem.rs b/src/models/algebraic/closest_vector_problem.rs index f707355a4..d306299f5 100644 --- a/src/models/algebraic/closest_vector_problem.rs +++ b/src/models/algebraic/closest_vector_problem.rs @@ -98,6 +98,42 @@ impl VarBounds { _ => None, } } + + /// Returns an exact bounded binary basis for offsets in this range. + /// + /// For a bounded variable with offsets `0..=hi-lo`, the returned weights + /// ensure that every bit-pattern reconstructs an in-range offset. Low-order + /// weights use powers of two; the final weight is capped so the maximum + /// reachable offset is exactly `hi-lo`. + pub(crate) fn exact_encoding_weights(&self) -> Vec { + let Some(num_values) = self.num_values() else { + panic!("CVP QUBO encoding requires finite variable bounds"); + }; + if num_values <= 1 { + return Vec::new(); + } + + let max_offset = (num_values - 1) as i64; + let num_bits = (usize::BITS - (num_values - 1).leading_zeros()) as usize; + let mut weights = Vec::with_capacity(num_bits); + + for bit in 0..num_bits.saturating_sub(1) { + weights.push(1_i64 << bit); + } + + let covered_by_lower_bits = if num_bits <= 1 { + 0 + } else { + (1_i64 << (num_bits - 1)) - 1 + }; + weights.push(max_offset - covered_by_lower_bits); + weights + } + + /// Returns the number of encoding bits needed for the exact bounded basis. + pub(crate) fn num_encoding_bits(&self) -> usize { + self.exact_encoding_weights().len() + } } /// Closest Vector Problem (CVP). @@ -175,6 +211,11 @@ impl ClosestVectorProblem { &self.bounds } + /// Returns the total number of bounded-encoding bits used by the QUBO form. + pub fn num_encoding_bits(&self) -> usize { + self.bounds.iter().map(VarBounds::num_encoding_bits).sum() + } + /// Convert a configuration (offsets from lower bounds) to integer values. fn config_to_values(&self, config: &[usize]) -> Vec { config diff --git a/src/rules/closestvectorproblem_qubo.rs b/src/rules/closestvectorproblem_qubo.rs new file mode 100644 index 000000000..69a2941b5 --- /dev/null +++ b/src/rules/closestvectorproblem_qubo.rs @@ -0,0 +1,197 @@ +//! Reduction from ClosestVectorProblem to QUBO. +//! +//! Encodes each bounded CVP coefficient with an exact in-range binary basis and +//! expands the squared-distance objective into a QUBO over those bits. + +#[cfg(feature = "example-db")] +use crate::export::SolutionPair; +use crate::models::algebraic::{ClosestVectorProblem, QUBO}; +use crate::reduction; +use crate::rules::traits::{ReduceTo, ReductionResult}; + +#[derive(Debug, Clone)] +struct EncodingSpan { + start: usize, + weights: Vec, +} + +/// Result of reducing a bounded ClosestVectorProblem instance to QUBO. +#[derive(Debug, Clone)] +pub struct ReductionCVPToQUBO { + target: QUBO, + encodings: Vec, +} + +impl ReductionResult for ReductionCVPToQUBO { + type Source = ClosestVectorProblem; + type Target = QUBO; + + fn target_problem(&self) -> &Self::Target { + &self.target + } + + /// Reconstruct the source configuration offsets from the encoded QUBO bits. + fn extract_solution(&self, target_solution: &[usize]) -> Vec { + self.encodings + .iter() + .map(|encoding| { + encoding + .weights + .iter() + .enumerate() + .map(|(offset, weight)| { + target_solution + .get(encoding.start + offset) + .copied() + .unwrap_or(0) + * weight + }) + .sum() + }) + .collect() + } +} + +#[cfg(feature = "example-db")] +fn canonical_cvp_instance() -> ClosestVectorProblem { + ClosestVectorProblem::new( + vec![vec![2, 0], vec![1, 2]], + vec![2.8, 1.5], + vec![ + crate::models::algebraic::VarBounds::bounded(-2, 4), + crate::models::algebraic::VarBounds::bounded(-2, 4), + ], + ) +} + +fn encoding_spans(problem: &ClosestVectorProblem) -> Vec { + let mut start = 0usize; + let mut spans = Vec::with_capacity(problem.num_basis_vectors()); + for bounds in problem.bounds() { + let weights = bounds + .exact_encoding_weights() + .into_iter() + .map(|weight| usize::try_from(weight).expect("encoding weights must be nonnegative")) + .collect::>(); + spans.push(EncodingSpan { start, weights }); + start += spans.last().expect("just pushed").weights.len(); + } + spans +} + +fn gram_matrix(problem: &ClosestVectorProblem) -> Vec> { + let basis = problem.basis(); + let n = basis.len(); + let mut gram = vec![vec![0.0; n]; n]; + for i in 0..n { + for j in i..n { + let dot = basis[i] + .iter() + .zip(&basis[j]) + .map(|(&lhs, &rhs)| lhs as f64 * rhs as f64) + .sum::(); + gram[i][j] = dot; + gram[j][i] = dot; + } + } + gram +} + +fn at_times_target(problem: &ClosestVectorProblem) -> Vec { + problem + .basis() + .iter() + .map(|column| { + column + .iter() + .zip(problem.target()) + .map(|(&entry, &target)| entry as f64 * target) + .sum() + }) + .collect() +} + +#[reduction(overhead = { num_vars = "num_encoding_bits" })] +impl ReduceTo> for ClosestVectorProblem { + type Result = ReductionCVPToQUBO; + + fn reduce_to(&self) -> Self::Result { + let encodings = encoding_spans(self); + let total_bits = encodings + .last() + .map(|encoding| encoding.start + encoding.weights.len()) + .unwrap_or(0); + let mut matrix = vec![vec![0.0; total_bits]; total_bits]; + + if total_bits == 0 { + return ReductionCVPToQUBO { + target: QUBO::from_matrix(matrix), + encodings, + }; + } + + let gram = gram_matrix(self); + let h = at_times_target(self); + let lowers = self + .bounds() + .iter() + .map(|bounds| { + bounds + .lower + .expect("CVP QUBO reduction requires finite lower bounds") + }) + .map(|lower| lower as f64) + .collect::>(); + let g_lo_minus_h = (0..self.num_basis_vectors()) + .map(|i| { + (0..self.num_basis_vectors()) + .map(|j| gram[i][j] * lowers[j]) + .sum::() + - h[i] + }) + .collect::>(); + + let mut bit_terms = Vec::with_capacity(total_bits); + for (var_index, encoding) in encodings.iter().enumerate() { + for &weight in &encoding.weights { + bit_terms.push((var_index, weight as f64)); + } + } + + for u in 0..total_bits { + let (var_u, weight_u) = bit_terms[u]; + matrix[u][u] = + gram[var_u][var_u] * weight_u * weight_u + 2.0 * weight_u * g_lo_minus_h[var_u]; + + for v in (u + 1)..total_bits { + let (var_v, weight_v) = bit_terms[v]; + matrix[u][v] = 2.0 * gram[var_u][var_v] * weight_u * weight_v; + } + } + + ReductionCVPToQUBO { + target: QUBO::from_matrix(matrix), + encodings, + } + } +} + +#[cfg(feature = "example-db")] +pub(crate) fn canonical_rule_example_specs() -> Vec { + vec![crate::example_db::specs::RuleExampleSpec { + id: "closestvectorproblem_to_qubo", + build: || { + crate::example_db::specs::rule_example_with_witness::<_, QUBO>( + canonical_cvp_instance(), + SolutionPair { + source_config: vec![3, 3], + target_config: vec![0, 0, 1, 0, 0, 1], + }, + ) + }, + }] +} + +#[cfg(test)] +#[path = "../unit_tests/rules/closestvectorproblem_qubo.rs"] +mod tests; diff --git a/src/rules/mod.rs b/src/rules/mod.rs index 68b87be76..a6b7b5fa1 100644 --- a/src/rules/mod.rs +++ b/src/rules/mod.rs @@ -7,6 +7,7 @@ pub use cost::{CustomCost, Minimize, MinimizeSteps, PathCostFn}; pub use registry::{ReductionEntry, ReductionOverhead}; pub(crate) mod circuit_spinglass; +mod closestvectorproblem_qubo; pub(crate) mod coloring_qubo; pub(crate) mod factoring_circuit; mod graph; @@ -89,6 +90,7 @@ pub use traits::{ReduceTo, ReductionAutoCast, ReductionResult}; pub(crate) fn canonical_rule_example_specs() -> Vec { let mut specs = Vec::new(); specs.extend(circuit_spinglass::canonical_rule_example_specs()); + specs.extend(closestvectorproblem_qubo::canonical_rule_example_specs()); specs.extend(coloring_qubo::canonical_rule_example_specs()); specs.extend(factoring_circuit::canonical_rule_example_specs()); specs.extend(knapsack_qubo::canonical_rule_example_specs()); diff --git a/src/unit_tests/models/algebraic/closest_vector_problem.rs b/src/unit_tests/models/algebraic/closest_vector_problem.rs index 7bb465789..d4b2cd5ff 100644 --- a/src/unit_tests/models/algebraic/closest_vector_problem.rs +++ b/src/unit_tests/models/algebraic/closest_vector_problem.rs @@ -56,6 +56,39 @@ fn test_cvp_dims() { assert_eq!(cvp.dims(), vec![5, 6]); // (-1..3)=5 values, (0..5)=6 values } +#[test] +fn test_cvp_num_encoding_bits_uses_var_bound_ranges() { + let basis = vec![vec![1, 0, 0], vec![0, 1, 0], vec![0, 0, 1]]; + let target = vec![0.0, 0.0, 0.0]; + let bounds = vec![ + VarBounds::bounded(-2, 4), + VarBounds::bounded(0, 5), + VarBounds::bounded(3, 3), + ]; + let cvp = ClosestVectorProblem::new(basis, target, bounds); + + assert_eq!(cvp.num_encoding_bits(), 6); +} + +#[test] +fn test_var_bounds_exact_encoding_weights_cap_non_power_of_two_ranges() { + let weights = VarBounds::bounded(0, 5).exact_encoding_weights(); + let represented_offsets: std::collections::BTreeSet = (0..(1usize << weights.len())) + .map(|mask| { + weights + .iter() + .enumerate() + .filter(|(bit, _)| ((mask >> bit) & 1) == 1) + .map(|(_, &weight)| weight) + .sum() + }) + .collect(); + + assert_eq!(weights, vec![1, 2, 2]); + assert_eq!(represented_offsets, (0..=5).collect()); + assert!(VarBounds::bounded(4, 4).exact_encoding_weights().is_empty()); +} + #[test] fn test_cvp_brute_force() { // b1=(2,0,0), b2=(1,2,0), b3=(0,1,2), target=(3,3,3) diff --git a/src/unit_tests/rules/closestvectorproblem_qubo.rs b/src/unit_tests/rules/closestvectorproblem_qubo.rs new file mode 100644 index 000000000..c1209d227 --- /dev/null +++ b/src/unit_tests/rules/closestvectorproblem_qubo.rs @@ -0,0 +1,85 @@ +use super::*; +use crate::models::algebraic::VarBounds; +use crate::rules::test_helpers::assert_optimization_round_trip_from_optimization_target; +use crate::solvers::BruteForce; + +fn canonical_cvp() -> ClosestVectorProblem { + ClosestVectorProblem::new( + vec![vec![2, 0], vec![1, 2]], + vec![2.8, 1.5], + vec![VarBounds::bounded(-2, 4), VarBounds::bounded(-2, 4)], + ) +} + +fn assert_close(actual: f64, expected: f64) { + assert!( + (actual - expected).abs() < 1e-9, + "expected {expected}, got {actual}" + ); +} + +#[test] +fn test_closestvectorproblem_to_qubo_closed_loop() { + let source = canonical_cvp(); + let reduction = ReduceTo::>::reduce_to(&source); + + assert_eq!(reduction.target_problem().num_vars(), 6); + assert_optimization_round_trip_from_optimization_target( + &source, + &reduction, + "ClosestVectorProblem->QUBO closed loop", + ); +} + +#[test] +fn test_closestvectorproblem_to_qubo_example_matrix_coefficients() { + let source = canonical_cvp(); + let reduction = ReduceTo::>::reduce_to(&source); + let qubo = reduction.target_problem(); + + assert_eq!(qubo.num_vars(), 6); + assert_close(*qubo.get(0, 0).expect("Q[0,0]"), -31.2); + assert_close(*qubo.get(0, 1).expect("Q[0,1]"), 16.0); + assert_close(*qubo.get(0, 2).expect("Q[0,2]"), 24.0); + assert_close(*qubo.get(1, 2).expect("Q[1,2]"), 48.0); + assert_close(*qubo.get(2, 5).expect("Q[2,5]"), 36.0); + assert_close(*qubo.get(5, 5).expect("Q[5,5]"), -73.8); +} + +#[test] +fn test_extract_solution_ignores_duplicate_exact_range_encodings() { + let reduction = ReduceTo::>::reduce_to(&canonical_cvp()); + + assert_eq!(reduction.extract_solution(&[1, 1, 0, 1, 1, 0]), vec![3, 3]); + assert_eq!(reduction.extract_solution(&[0, 0, 1, 0, 0, 1]), vec![3, 3]); +} + +#[cfg(feature = "example-db")] +#[test] +fn test_closestvectorproblem_to_qubo_canonical_example_spec() { + let spec = canonical_rule_example_specs() + .into_iter() + .find(|spec| spec.id == "closestvectorproblem_to_qubo") + .expect("missing canonical ClosestVectorProblem -> QUBO example spec"); + let example = (spec.build)(); + + assert_eq!(example.source.problem, "ClosestVectorProblem"); + assert_eq!(example.target.problem, "QUBO"); + assert_eq!(example.target.instance["num_vars"], 6); + assert_eq!(example.solutions[0].source_config, vec![3, 3]); + assert_eq!(example.solutions[0].target_config, vec![0, 0, 1, 0, 0, 1]); +} + +#[test] +fn test_duplicate_target_encodings_have_equal_qubo_value() { + let reduction = ReduceTo::>::reduce_to(&canonical_cvp()); + let qubo = reduction.target_problem(); + let solver = BruteForce::new(); + let best = solver.find_all_best(qubo); + + assert!(best.contains(&vec![0, 0, 1, 0, 0, 1]) || best.contains(&vec![1, 1, 0, 1, 1, 0])); + assert_close( + qubo.evaluate(&[0, 0, 1, 0, 0, 1]), + qubo.evaluate(&[1, 1, 0, 1, 1, 0]), + ); +}