Fix #125: [Rule] SubsetSum to ClosestVectorProblem (#709)

* Add plan for #125: [Rule] SubsetSum to ClosestVectorProblem

* Implement #125: [Rule] SubsetSum to ClosestVectorProblem

* chore: remove plan file after implementation

---------

Co-authored-by: Xiwei Pan <90967972+isPANN@users.noreply.github.com>

Author: GiggleLiu

docs/paper/reductions.typ

@@ -4586,6 +4586,50 @@ where $P$ is a penalty weight large enough that any constraint violation costs m
   _Solution extraction._ For each $i$: if $y_i$ is selected ($x_(2i) = 1$), set $x_i = 1$; if $z_i$ is selected ($x_(2i+1) = 1$), set $x_i = 0$.
 ]
 
+#{
+  let ss-cvp = load-example("SubsetSum", "ClosestVectorProblem")
+  let ss-cvp-sol = ss-cvp.solutions.at(0)
+  let ss-cvp-sizes = ss-cvp.source.instance.sizes
+  let ss-cvp-target = ss-cvp.source.instance.target
+  let ss-cvp-basis = ss-cvp.target.instance.basis
+  let ss-cvp-target-vec = ss-cvp.target.instance.target
+  let ss-cvp-n = ss-cvp-sizes.len()
+  let ss-cvp-x = ss-cvp-sol.target_config
+  let to-mat(m) = math.mat(..m.map(row => row.map(v => $#v$)))
+  [
+    #reduction-rule("SubsetSum", "ClosestVectorProblem",
+      example: true,
+      example-caption: [#ss-cvp-n elements, target sum $B = #ss-cvp-target$],
+      extra: [
+        *Step 1 -- Source instance.* The canonical Subset Sum instance has sizes $(#ss-cvp-sizes.map(str).join(", "))$ and target $B = #ss-cvp-target$.
+
+        *Step 2 -- Build the lattice.* The reduction creates the basis
+        $ bold(B) = #to-mat(ss-cvp-basis) $
+        together with target $ bold(t) = (#ss-cvp-target-vec.map(str).join(", "))^top $
+        and binary bounds $x_i in {0,1}$ for all $#ss-cvp-n$ coordinates.
+
+        *Step 3 -- Verify the canonical witness.* The fixture stores $bold(x) = (#ss-cvp-x.map(str).join(", "))$, which selects sizes $3$ and $8$ and therefore satisfies $3 + 8 = #ss-cvp-target$. Since $bold(B) bold(x) = (1, 0, 0, 1, #ss-cvp-target)^top$, the difference vector is $(0.5, -0.5, -0.5, 0.5, 0)^top$ and the Euclidean distance is $sqrt(#ss-cvp-n / 4) = 1$.
+
+        *Witness semantics.* The example DB stores one canonical minimizer. This source instance also has another satisfying subset, $(1, 1, 1, 0)$, so the reduction has multiple optimal CVP witnesses even though only one is serialized.
+      ],
+    )[
+      Classical lattice embedding for Subset Sum following Lagarias and Odlyzko @lagarias1985, with the $1/2$-target CVP formulation in the style of Coster et al. @coster1992. For an instance with $n$ elements, the reduction produces $n$ basis vectors in ambient dimension $n + 1$: the first $n$ coordinates enforce binary structure and the last coordinate records the subset sum error.
+    ][
+      _Construction._ Given sizes $s_0, dots, s_(n-1) in ZZ^+$ and target $B in ZZ^+$, define one basis vector per element:
+      $ bold(b)_i = bold(e)_i + s_i bold(e)_(n+1) $
+      for $i in {0, dots, n-1}$. Equivalently, the basis matrix has columns $bold(b)_0, dots, bold(b)_(n-1)$, so its first $n$ rows form the identity matrix and its last row is $(s_0, dots, s_(n-1))$. Set the target vector to
+      $ bold(t) = (1/2, dots, 1/2, B)^top $
+      and restrict every CVP variable to $x_i in {0, 1}$.
+
+      _Correctness._ ($arrow.r.double$) If $bold(x) in {0,1}^n$ is a satisfying Subset Sum solution, then $sum_i s_i x_i = B$ and
+      $ norm(bold(B) bold(x) - bold(t))_2^2 = sum_(i=0)^(n-1) (x_i - 1/2)^2 + (sum_i s_i x_i - B)^2 = n/4. $
+      Hence every satisfying subset becomes a CVP solution at distance $sqrt(n / 4)$. ($arrow.l.double$) Conversely, binary bounds force every CVP candidate to lie in ${0,1}^n$. The first $n$ coordinates always contribute exactly $n/4$ to the squared distance, so a CVP minimizer attains distance $sqrt(n/4)$ if and only if the last coordinate contributes $0$, i.e. $sum_i s_i x_i = B$. When the Subset Sum instance is unsatisfiable, every binary vector has strictly larger distance.
+
+      _Solution extraction._ Return the binary CVP vector unchanged.
+    ]
+  ]
+}
+
 #reduction-rule("ILP", "QUBO")[
   A binary ILP optimizes a linear objective over binary variables subject to linear constraints. The penalty method converts each equality constraint $bold(a)_k^top bold(x) = b_k$ into the quadratic penalty $(bold(a)_k^top bold(x) - b_k)^2$, which is zero if and only if the constraint is satisfied. Inequality constraints are first converted to equalities using binary slack variables with powers-of-two coefficients. The resulting unconstrained quadratic over binary variables is a QUBO whose matrix $Q$ combines the negated objective (as diagonal terms) with the expanded constraint penalties (as a Gram matrix $A^top A$).
 ][

docs/paper/references.bib

@@ -143,6 +143,28 @@ @inproceedings{karp1972
   pages     = {85--103}
 }
 
+@article{lagarias1985,
+  author  = {Jeffrey C. Lagarias and Andrew M. Odlyzko},
+  title   = {Solving Low-Density Subset Sum Problems},
+  journal = {Journal of the ACM},
+  volume  = {32},
+  number  = {1},
+  pages   = {229--246},
+  year    = {1985},
+  doi     = {10.1145/2455.2461}
+}
+
+@article{coster1992,
+  author  = {Matthijs J. Coster and Antoine Joux and Brian A. LaMacchia and Andrew M. Odlyzko and Claus-Peter Schnorr and Jacques Stern},
+  title   = {Improved Low-Density Subset Sum Algorithms},
+  journal = {Computational Complexity},
+  volume  = {2},
+  number  = {2},
+  pages   = {111--128},
+  year    = {1992},
+  doi     = {10.1007/BF01201999}
+}
+
 @inproceedings{cook1971,
   author    = {Stephen A. Cook},
   title     = {The Complexity of Theorem-Proving Procedures},

src/rules/mod.rs

@@ -37,6 +37,7 @@ pub(crate) mod sat_minimumdominatingset;
 mod spinglass_casts;
 pub(crate) mod spinglass_maxcut;
 pub(crate) mod spinglass_qubo;
+pub(crate) mod subsetsum_closestvectorproblem;
 #[cfg(test)]
 pub(crate) mod test_helpers;
 mod traits;
@@ -115,6 +116,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
     specs.extend(sat_minimumdominatingset::canonical_rule_example_specs());
     specs.extend(spinglass_maxcut::canonical_rule_example_specs());
     specs.extend(spinglass_qubo::canonical_rule_example_specs());
+    specs.extend(subsetsum_closestvectorproblem::canonical_rule_example_specs());
     specs.extend(travelingsalesman_qubo::canonical_rule_example_specs());
     #[cfg(feature = "ilp-solver")]
     {

src/rules/subsetsum_closestvectorproblem.rs

@@ -0,0 +1,83 @@
+//! Reduction from Subset Sum to Closest Vector Problem.
+
+use crate::models::algebraic::{ClosestVectorProblem, VarBounds};
+use crate::models::misc::SubsetSum;
+use crate::reduction;
+use crate::rules::traits::{ReduceTo, ReductionResult};
+use num_bigint::BigUint;
+use num_traits::ToPrimitive;
+
+/// Result of reducing SubsetSum to ClosestVectorProblem.
+#[derive(Debug, Clone)]
+pub struct ReductionSubsetSumToClosestVectorProblem {
+    target: ClosestVectorProblem<i32>,
+}
+
+impl ReductionResult for ReductionSubsetSumToClosestVectorProblem {
+    type Source = SubsetSum;
+    type Target = ClosestVectorProblem<i32>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        target_solution.to_vec()
+    }
+}
+
+fn biguint_to_i32(value: &BigUint) -> i32 {
+    value
+        .to_i32()
+        .expect("SubsetSum -> ClosestVectorProblem requires all sizes and target to fit in i32")
+}
+
+#[reduction(
+    overhead = {
+        ambient_dimension = "num_elements + 1",
+        num_basis_vectors = "num_elements",
+    }
+)]
+impl ReduceTo<ClosestVectorProblem<i32>> for SubsetSum {
+    type Result = ReductionSubsetSumToClosestVectorProblem;
+
+    fn reduce_to(&self) -> Self::Result {
+        let n = self.num_elements();
+        let mut basis = Vec::with_capacity(n);
+        for (i, size) in self.sizes().iter().enumerate() {
+            let mut column = vec![0i32; n + 1];
+            column[i] = 1;
+            column[n] = biguint_to_i32(size);
+            basis.push(column);
+        }
+
+        let mut target = vec![0.5; n];
+        target.push(biguint_to_i32(self.target()) as f64);
+
+        ReductionSubsetSumToClosestVectorProblem {
+            target: ClosestVectorProblem::new(basis, target, vec![VarBounds::binary(); n]),
+        }
+    }
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    use crate::export::SolutionPair;
+
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "subsetsum_to_closestvectorproblem",
+        build: || {
+            crate::example_db::specs::rule_example_with_witness::<_, ClosestVectorProblem<i32>>(
+                SubsetSum::new(vec![3u32, 7, 1, 8], 11u32),
+                SolutionPair {
+                    source_config: vec![1, 0, 0, 1],
+                    target_config: vec![1, 0, 0, 1],
+                },
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/subsetsum_closestvectorproblem.rs"]
+mod tests;

src/unit_tests/rules/subsetsum_closestvectorproblem.rs

@@ -0,0 +1,79 @@
+use super::*;
+use crate::models::algebraic::{ClosestVectorProblem, VarBounds};
+use crate::rules::test_helpers::assert_satisfaction_round_trip_from_optimization_target;
+use crate::solvers::{BruteForce, Solver};
+use crate::traits::Problem;
+use crate::types::SolutionSize;
+use std::collections::HashSet;
+
+#[test]
+fn test_subsetsum_to_closestvectorproblem_closed_loop() {
+    let source = SubsetSum::new(vec![3u32, 7, 1, 8], 11u32);
+    let reduction = ReduceTo::<ClosestVectorProblem<i32>>::reduce_to(&source);
+    let target = reduction.target_problem();
+
+    assert_eq!(target.num_basis_vectors(), 4);
+    assert_eq!(target.ambient_dimension(), 5);
+    assert_eq!(target.bounds(), &[VarBounds::binary(); 4]);
+
+    assert_satisfaction_round_trip_from_optimization_target(
+        &source,
+        &reduction,
+        "SubsetSum -> ClosestVectorProblem closed loop",
+    );
+}
+
+#[test]
+fn test_subsetsum_to_closestvectorproblem_structure() {
+    let source = SubsetSum::new(vec![3u32, 7, 1, 8], 11u32);
+    let reduction = ReduceTo::<ClosestVectorProblem<i32>>::reduce_to(&source);
+    let target = reduction.target_problem();
+
+    assert_eq!(target.basis()[0], vec![1, 0, 0, 0, 3]);
+    assert_eq!(target.basis()[1], vec![0, 1, 0, 0, 7]);
+    assert_eq!(target.basis()[2], vec![0, 0, 1, 0, 1]);
+    assert_eq!(target.basis()[3], vec![0, 0, 0, 1, 8]);
+    assert_eq!(target.target(), &[0.5, 0.5, 0.5, 0.5, 11.0]);
+}
+
+#[test]
+fn test_subsetsum_to_closestvectorproblem_issue_example_minimizers() {
+    let source = SubsetSum::new(vec![3u32, 7, 1, 8], 11u32);
+    let reduction = ReduceTo::<ClosestVectorProblem<i32>>::reduce_to(&source);
+    let target = reduction.target_problem();
+    let solutions: HashSet<Vec<usize>> = BruteForce::new()
+        .find_all_best(target)
+        .into_iter()
+        .collect();
+
+    let expected: HashSet<Vec<usize>> = [vec![1, 0, 0, 1], vec![1, 1, 1, 0]].into_iter().collect();
+    assert_eq!(solutions, expected);
+
+    for solution in &solutions {
+        assert_eq!(target.evaluate(solution), SolutionSize::Valid(1.0));
+    }
+}
+
+#[test]
+fn test_subsetsum_to_closestvectorproblem_unsatisfiable_instance() {
+    let source = SubsetSum::new(vec![2u32, 4, 6], 5u32);
+    let reduction = ReduceTo::<ClosestVectorProblem<i32>>::reduce_to(&source);
+    let target = reduction.target_problem();
+    let best = BruteForce::new()
+        .find_best(target)
+        .expect("unsatisfiable instance should still have a best CVP assignment");
+
+    match target.evaluate(&best) {
+        SolutionSize::Valid(value) => assert!(value > (source.num_elements() as f64).sqrt() / 2.0),
+        SolutionSize::Invalid => panic!("CVP solution should be valid"),
+    }
+}
+
+#[test]
+#[should_panic(
+    expected = "SubsetSum -> ClosestVectorProblem requires all sizes and target to fit in i32"
+)]
+fn test_subsetsum_to_closestvectorproblem_panics_on_large_coefficients() {
+    let source = SubsetSum::new(vec![(i32::MAX as u64) + 1], 1u64);
+    let _ = ReduceTo::<ClosestVectorProblem<i32>>::reduce_to(&source);
+}